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ABSTRACT 


A theory is developed for a digital equalizer for 
use in reducing intersymbol interference (ISI) on high- 
speed data communi cations channels. The equalizer is 
initialized with a single isolated transmitter pulse, ' 
provided the signal-to-noise ratio (SNR) is not unusually 
low, then switches to a decision-directed, on-line mode 
of operation that allows tracking of channel variations. 
Conditions for optimal tap-gain settings are obtained first 
for a transversal equalizer structure by using a mean- 
squared error (MSE) criterion, a first-order gradient 
algorithm to determine the adjustable equalizer tap-gains, 
and a sequence of isolated initializing pulses. Since the 
rate of tap-gain convergence depends on the eigenvalues of 
a channel-output correlation matrix A convergence can be 
improved by making a linear transformation on A to obtain 
a new correlation matrix A. The corresponding change in 
equalizer structure results in a set of parallel trans- 
versal filter-sections, the weighted outputs of which are 
summed to produce the equalizer output. It is shown that 
by properly selecting the tap-gains of each transversal 
filter-section, the equalizer can be adjusted to minimize 
the MSE on the first iteration of the algorithm. A method 
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is proposed for making this selection using the pseudoinverse 
of A to obtain a "best" approximation to an orthogonality 
condition. By using a Hadamard matrix it is shown that A 
is similar to a projection matrix. Thus all eigenvalues 
are either 0 to 1, which is the desired condition for 
initialization in one algorithm iteration. 

When the channel SNR is low initialization can be 
improved by averaging several initializing pulses in order 
to better identify the channel. Next a simple approxi- 
mation to the error criterion is used to obtain an on-line 
algorithm. for noisy channels; it makes small adjustments 
to the filter-section output weights to compensate for 
channel noise and time-variations, as well as poor initiali- 
zation. Convergence of this on-line algorithm is shown in 
a statistical sense under a stated hypothesis, and 
suggestions are made for detection and correction of poor 
channel tracking. Finally the equalizer design is evaluated 
by means of a digital computer simulation and the results 
are presented and discussed for three baseband channel 
models. The interaction among several key equalizer para- 
meters is studied for a wide range of SNR. 
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CHAPTER 1 


INTRODUCTION 

Motivated by the need to transmit binary data 
directly between computers, there has been a steady increase 
in the use of digital techniques for data transmission. 

This is due in large part to the profound advances in solid- 
state technology which have resulted in small, reliable, 
fast, and inexpensive digital components. Attempts to 
increase the transmission rate of digital data over band- 
limited channels are often limited by interference among 
successive data symbols, an effect known as intersymbol 
interference (ISI) . The reduction of ISI by signal 
processing techniques is known as equalization . 

1.1. THE EQUALIZATION OF DIGITAL CHANNELS 

The basic principles ;.of digital communications are 
well documented in standard texts [1]. Such systems always 
contain unwanted signals which are usually classified 
collectively as noise . For this study it is convenient to 
divide noise into two subcategories? background noise caused 
by external sources such as thermal excitation and ISI 
caused by the overlap in time of the received signals. As 
the data transmission rate is increased for a fixed channel 
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bandwidth, the ISI becomes increasingly important for the 
following reason. Over the nominal channel passband all 
practical channels have nonlinear phase vs. frequency 
characteristics which cause different frequency components 
of a transmitted signal to be delayed by different amounts. 
This results in a transmitted pulse being spread out, or 
dispersed , in time. This time dispersion is further 
aggravated by the fact that all physical channels are 
bandlimited and thus attenuate the high frequency components 
of transmitted signals. At the channel output these 
successive time dispersed signals overlap to give a composite 
signal equal to their algebraic sum. Fig. 1.1 shows such a 
situation, with the circled points indicating those values 
measured by a receiver that takes a single sample per bit. 

In the figure time has been normalized, and T^ represents 
the normalized time delay of the transmitted pulse caused 
by the channel. The channel output sample at t = + 2, 

for example, is closer to unity than zero? therefore, the 
zero transmitted at t = 2 would be mistaken for a one, 
resulting in a message error. Errors can also be caused 
by background noise or time variation of channel parameters. 
Some ways to reduce the error rate include raising the 
signal-to-noise ratio (SNR) , introducing redundancy through 
coding, and the use of pulse-shaping filters to reshape the 
transmitted or received signals. This latter approach is 
called channel equalization, and the filter is termed an 
equalizer . 
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A number of classification schemes have been proposed 
for equalizers; the most important for this study being 
analog or digital according to the type of signal being 
processed. Digital equalizers may be further classified by 
the duration of their response to a standard test signal 
called the unit pulse , i.e. , a binary one preceded and 
followed by strings of binary zeros. The unit pulse is 
the discrete analog of the unit impulse for continuous 
systems, and is in fact often referred to as the "discrete 
unit impulse." If the theoretical unit-pulse response goes 
to zero after a finite time and remains there, then the 
equalizer has a finite impulse response (FIR) , otherwise it 
has an infinite impulse response (HR) . Other classification 
schemes are based on linearity (linear, nonlinear) , feedback 
(recursive, nonrecursive), signal flow (parallel, series, 
cascade) , stationarity of characteristics (time-varying, 
time -invariant) , and historical nomenclature tfrequency- 
sampling, transversal, sampled data) . A thorough discussion 
of these classifications is given by Rabiner, et al. [2]. 

Of particular interest here is the transversal filter 
structure , shown in Fig. 1.2. The blocks represent delay 
elements of duration T s, such as digital shift registers 
or analog tapped delay line stages, and the numbers {c n > 
are called tap-gains . It will be convenient to mathemati- 
cally represent the T s delay by the z-transform delay 

_ i j. 

operator z which is discussed in A. 4. Transversal 
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The notation A. 4 refers to Section 4 of App. A. 
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equalizers have a number of nice properties such as good 
stability characteristics, ease of implementation, linearity, 
and mathematical tractability . Consequently, they have 
proven particularly useful in equalizing ISI. 

If an equalizer is set at some initial time, e.g., 
when installed, and not changed thereafter, then it is 
called a preset or fixed equalizer. Fixed equalizers often 
prove inadequate in modern digital communications networks 
having switching between many branches, atmospheric trans- 
mission links, wide temperature fluctuations, and component 
aging. Such applications require an equalizer capable of 
being adjusted either during normal message transmission or 
at selected breaks between messages. Although the term 
adaptive equalization has sometimes been used for both of 
these cases, it will be reserved herein to describe only 
the former, while the latter case will be called automatic 
equalization [3] . Both adaptive and automatic equalizers 
will be referred to as self-adj usting . 

Since a separate feedback channel from receiver to 
transmitter is usually infeasible in complex communication 
systems, some alternative method is necessary to produce an 
error signal for use by the equalizer. In adaptive 
operation this can be done by employing a dec is ion -directed 
scheme. It is assumed that a high percentage of the output 
decisions are correct so that the sequence of output 
decisions closely approximates the transmitted sequence. 
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The error signal is then defined as the difference between 

/ 

signals before and after the decisions are made. The 
limitation of this scheme is obvious? if the equalizer is 
poorly adjusted initially, then the output sequence contains 
so many errors that the error signal becomes meaningless. 
Consequently, it is usually necessary to transmit an initial 
known training or test sequence which can be reproduced at 
the receiver. Comparison with the equalizer output before 
the decision device then gives the desired error signal. 

This known test sequence, transmitted during the initialization 
or training period is frequently selected to be a string of 
unit pulses with sufficient spacing between pulses to prevent 
the received waveforms from overlapping. After each such 
isolated test pulse an incremental adjustment of a particular 
set of equalizer parameters is made according to an iterative 
procedure called an algorithm . Assuming that the algorithm 
is convergent, the training period is continued until 
initialization is achieved, at which time the message 
transmission commences. In applications such as networks of 
sensors or polling systems where the messages are short, 
this training period can be considerably longer than the 
message length. Furthermore, many common equalizer 
structures require that the adjustment increments be small 
to insure convergence of the algorithm, thereby necessi- 
tating lengthy test sequences. 
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The problem to be studied herein is the development 
of an equalizer with a short initialization period that can 
be operated in an adaptive manner after initialization. 

It should have practical application to noisy time-varying 
networks that include switching. It is assumed throughout 
that sufficient digital computing equipment is available, 
such as a hardwired logic network, a minicomputer, or a 
time-shared general purpose computing system. 

1.2 BACKGROUND AND ORGANIZATION 

Some familiarity with the theory of modern digital 
communications in both the time and frequency domains is 
assumed herein, and a basic knowledge of applied mathe- 
matics is essential. The discrete nature of digital signals 
makes the z-transform a valuable tool, and some essential 
properties are given in App . A. This appendix also contains 
some fundamentals of linear operators which are then 
applied to the treatment of selected topics from the theory 
of vectors and matrices. Familiarity with probability and 
statistics is also desirable. A number of references are 
cited in the appendices both for completeness and as a 
source of proofs omitted from App. B. Finally, a reading 
knowledge of the FORTRAN programming language will be 
helpful in understanding the computer simulation listings 
found in App. C. 
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In addition to some introductory material on digital 
equalizers, this first chapter also contains a review of 
related work found in the literature. In Chapter 2 the 
initialization theory is developed for an automatic digital 
equalizer under the assumption that there is no background 
noise, i.e., the channel is deterministic. The equalizer 
uses a first-order gradient algorithm and a short test 
sequence to achieve initialization on the first iteration of 
the algorithm. A number of proofs relating to this chapter 
were relegated to App. B in order to preserve continuity. 
Chapter 3 first considers the effect of noise on the 
previously developed initialization theory. Then an 
approximation to the gradient algorithm is introduced which 
allows the equalizer to make on-line adjustments for time- 
variation of the channel parameters and noise statistics. 

The resulting adaptive equalizer is evaluated by means of 
a digital computer simulation, with the modeling assumptions 
presented and discussed. A complete description of the 
computer simulation program including a listing is found in 
App. C. The concluding chapter. Chapter 4, includes some 
practical considerations for implementation of the equalizer 
design developed in the preceding chapters and itemizes some 
additional areas for further investigation. 
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1.3 SURVEY OF RELATED LITERATURE 

Equalization has been used on analog communication 
channels for many years to improve the channel amplitude 
vs. frequency characteristics. Such equalizers were 
generally of the fixed type, and the resulting channel phase 
vs. frequency characteristics were usually unimportant since 
the human ear is not very sensitive to the resulting delay 
distortion. However, with the advent in the 1950 ' s of 
extensive digital traffic over the increasingly complex 
existing networks, it became more and more important to 
reduce ISI by making the phase vs. frequency characteristic 
as linear as possible in the channel passband. These 
attempts at increasing the data rates by decreasing ISI 
stimulated considerable interest in the use of transversal 
equalizers employing tapped delay lines and analog multipliers , 
e.g., attenuators, for the tap-gain multipliers. The first 
noteworthy attempt to make such equalizers self-adjusting 
was by Lucky [4] in 1965. By defining the peak distortion 
error criterion as a measure of ISI he obtained conditions 
for the optimal tap-gain settings using a series of isolated 
initializing pulses prior to message transmission in con- 
junction with a steepest-descent gradient technique. 

By computing the approximate gradient with only the polarity 
of the variables involved he was able to realize the 
equalizer using only digital logic. Lucky subsequently 
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combined the initialization procedure with a decision- 
directed on-line adaptive scheme for slowly-varying channels, 
and was again able to realize the equalizer with digital 
circuits [3]. Although these equalizers performed, well when 
initialized, the initialization procedure generally required 
a large number of initializing pulses and did not converge at 
all when the peak distortion was greater than unity. Such 
equalizers were simple to implement and were inherently 
stable. They were named zero-forcing, since the unit-pulse 
response of the equalized channel was made zero at a number 
of points adjacent to the main pulse, with the number of 
such zeros dependent on the equalizer length. 

In 1966 Widrow [5] introduced a transversal equalizer 
based on optimization with the mean-squared error criterion 
[6], and made it adaptive by using a decision-directed 
scheme and an appropriate "noisy" gradient estimate to 
update the tap-gains. Although this equalizer required 
linear multipliers rather than the less complex logic for 
computing the gradient, the speed of convergence was much 
faster. Initialization still required a large number of 
isolated initializing pulses. The application of this 
equalizer was first considered by Lucky and Rudin 17] , and 
a detailed mathematical analysis of its convergence 
properties was done by Gersho £ 8 J . An interesting all- 
digital serial realization of this type of equalizer was 
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developed by Lender [9]. Hirsch and Wolf [103 described an 
actual equalizer developed for Bell Labs using integrated 
circuits and field-effect transistors for the tap-gain 
multipliers. General treatments of these linear self- 
adjusting equalizers have been done [11,123 , and the effects 
of quantization noise in digital implementations has also 
been studied [13] . 

In addition to the linear transversal equalizers, a 
number of other approaches have been made to the design 
and improvement of digital equalizers. One method 1143 used 
the recursive frequency-sampling structure with a special 
initializing sequence that allowed each two-pole filter to 
be' adjusted independently. Use of the Viterbi algorithm 115] 
for on-line channel identification for equalization was 
discussed by Magee and Prokias [16] , and the somewhat 
related technique of dynamic programming has been suggested 

[17] to improve the initialization speed of an equalizer using 
a first-order gradient algorithm. By using the technique 

of quasilinearization the Kalman filter theory has been 
applied to estimate both the channel and the channel input 

[18] . Another approach to digital equalization uses the 
feedback of previous decisions through a transversal filter 
to compensate for the postcursors of the current pulse. 

This method is called decision-feedback [19,20], and is 
particularly good where the channel distortion is large with 
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relatively few precursors. However, decision errors tend 
to propagate and cause error bursts. 

Throughout the development of digital equalizers 
there have been continuing attempts to reduce the time 
required for initialization. For short messages the initiali- 
zation times were often excessive, e.g., an equalizer using 
the mean-squared error criteria was reported 110] to require 
several thousand isolated initializing pulses, and a 
considerably larger number of pulses was required when the 
equalization was done using the more slowly convergent zero- 
forcing equalizer. A number of methods were proposed to 
speed up this convergence, including the use of second-order 
gradient algorithms [21] , use of different error criteria 
[22] , and the periodic adjustment of the convergence factor 
in first-order algorithms [23,24]. This latter approach 
included the so-called "gearshift" method which used large 
initial adjustments to approach the optimal tap-gain value, 
and then changed to small adjustment increments for final 
accuracy. In 1971 Chang [25] proposed a new equalizer 
structure based on a set of parallel transversal filters. 

He showed that when the channel amplitude vs. frequency 
characteristic was accurately known, initialization was 
possible in one iteration of the tap-gain adjustment 
algorithm, and also that fast initialization was possible 
for a common class of partial response systems when the 



14 


channel was estimated. A similar structure has also been 
used to study the case when there are a small number or 
possible channels, each with a priori known statistics, to 
be equalized with a single equalizer [26]. There is 
currently a continuing interest in equalizer design, and 
the steady improvements in speed and availability of solid- 
state digital devices indicates that more sophisticated 
designs will soon be practical. 

1 . 4 COMMENTS ON NOTATION 

The theory of digital filters and equalizers has 
developed from a number of mathematical and engineering 
disciplines, e.g., networks, communications, control 
systems , optimization, and mathematical programming. 
Consequently, the notation and terminology found in the 
literature vary considerably. The most significant attempt 
at terminology standardization to date is found in 
Terminology in Digital Signal Processing [2] compiled by 
the committee on Digital Signal Processing of the IEEE 
Group on Audio and Electroacoustics. 

A list of symbols used herein is found on page 
iv/ i ‘however, some additional comments 'will also prove 

useful. All vectors are taken to be column vectors, and 
vector and matrix variables are indicated by an underline. 

To avoid confusion brackets are used to index the algorithm 
iterations. Several types of sequence notation are also 



used. The sequences in Chapter 2 generally relate to the 
first iteration of the tap-gain adjustment algorithm and 
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are indexed by subscripts, e.g. {x^}, while those in 
Chapter 3 relate one-to-one with the algorithm iterations 
so the index is placed in brackets, e.g. {xlkj}. Sequences 
are also represented by the corresponding z-transform, 
e.g. X(z), in the manner discussed in A. 4. When integrating 
in the z-domain the unit circle is always taken as the con- 
tour of integration. When evaluating these contour integrals 
it is often useful to make the substitution z = e^* 0 , where 
w is the normalized frequency in radians per sample, i.e., 

0 ^ b> < 2 tt . The resulting function will be denoted in the 
conventional manner, e.g., X(joa) vice X(e J j. Other notation 
that is specific to a particular section is defined therein 
and is not discussed here. This especially true of the 
appendices . 
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A THEORY FOR ONE-STEP EQUALIZATION 
IN THE ABSENCE OF NOISE 

After defining the channel model and the equalizer 
structure, an algorithm for equalizer adjustment is 
selected. The characteristics that control convergence of 
the transversal equalizer are studied, and from them a 
transformation is developed which insures convergence on 
the first iteration of the algorithm provided that the 
channel is noise free. 

Before proceeding, however, it is necessary to 
briefly examine the overall digital communication system 
shown in Fig. 2.1. The transmitter produces a sequence of 
digital symbols {d^}, taken here to be binary. This 
sequence is sent through a digital channel which distorts 
the input to produce the channel output sequence {x^}. 

The digital channel is described by its discrete unit-pulse 
response {g m ), but will usually include an analog channel 
such as a cable or radio link. It is then necessary to 
use a modulator at the analog channel input to map the 
digital sequence into a train of short analog waveforms, 
and a demodulator at the channel output to perform the 
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Fig. 2.1 Digital Communication System 
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inverse operation. In order to compensate for the effect 
of the channel an equalizer is employed at the channel 
output, and finally a decision device is used to obtain 
an estimate of the transmitted sequence based on 
observation of the equalizer output sequence {y^}. Note 
that each block of Fig. 2.1 involves only digital inputs 
and outputs, hence these signals can be represented in the 
time domain by sequences , e. g. , {xj^}, and in the transform 
domain by the respective z-transform 

OD 

X<2> - E x m z " m • 

m=0 

It will be assumed that all necessary timing and 
synchronization information is available, and that the 
channel time delay is known. These problems are 
extensively covered in the literature [2 7 #2 8 ] ■ The time 
delay will be absorbed in the notation; for example, the 
received signal x^ will correspond to the transmitted 
signal although it actually occurs T^-s later. 

2.1 CHANNEL MODEL AND THE GENERALIZED EQUALIZER STRUCTURE 

The type of channel to be considered is a linear 
digital channel at baseband having a finite-length 
response {g Q , 9}. ' ’ * * ^G-l^ - to the unit pulse, where G is 
the sequence length. The channel will in general be time 
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varying due to switching and other background noise. For 
this chapter it will be assumed that the channel is free of 
stochastic noise, and that no switching occurs during the 
initialization period. Furthermore, other channel variations 
will be assumed slow in relation to the response time for a 
single pulse. These assumptions allow channel time 
variations to be neglected during the one-step initialization 
period. The effect of stochastic noise will be considered 
in Chapter 3. In addition to the unit-pulse response, the 
channel can be described by a transfer function 

G-l 

GCs) = £ g m z -m , (2-1) 

m=0 

where G(z) is the z -transform of {go} . 

m 

A fixed equalizer can in theory be located at any 
point of a linear channel. However, a self-adjusting 
equalizer is usually located at the channel output since 
it must obtain sufficient information about the channel 
to initialize itself and then maintain a suitable setting 
as the channel and noise characteristics change. When this 
channel identification is performed using a known test 
sequence or by a decision-directed scheme, the result is 
obtained at the channel output, and so any other equalizer 
placement necessitates a highly reliable feedback channel. 
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Many linear automatic equalizer designs can be put 
in the parallel form shown in Fig. 2.2, which is referred 
to as the generalized equalizer structure I7].. The blocks 
{Fft(z)} will be called the filter-sections since each 
transfer function F m (z) represents a digital filter. When 
the channel output X(z) is applied to the filter-sections, 
the equalizer output Y(z) is obtained by multiplying each 
filter-section output {V^(z) }by the corresponding real -valued 
tap-gain c m and then summing over the M weighted outputs 
{c m Vm(z)}. The equalizer transfer function H(z) = Y(z)/X(z) 
is determined as follows: 

M-l 

H(z) = Xlz) E C m V m (z) = ^ *< z >' 

m=0 

where e = (fig ,c x , . . .c M _ 1 _) T and F(z) = (Fq (z) ,F 1 (z) , 
...F M _^(z)) . For self-adjusting operation, the tap-gains 
are made adjustable. A device called an adapter seeks to 
minimize the ISI by iteratively computing a set of varying 
tap-gain increments {Ac^tk]} using a predetermined algorithm. 
The tap-gains are then adjusted at iteration according to 

c m [k+l] = c m [k] + 4c m [k] , xn=0,l, . . .M-l. 

2.2 THE ALGORITHM AND ERROR CRITERION 

The usefulness of the first-order gradient algorithm 
for self-adjusting equalizers is well established in the 
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literature I5,7,8J . It has the vector form 

c [k+1] = c[k] - 1/2 y 7J(c[k]), (2.2) 

where c[k] is the M-vector of tap-gains at iteration k and 
the corresponding tap-gain increment vector is given by 
Ac [kj = c([k+l] - c[k] . The real scalar y is called the 
convergence factor » It scales the tap-gain increments and 
thereby determines whether the sequence of tap-gain vectors 
{c|k] } converges to some final value, and if so, the rate 
of this convergence. The differentiable, real-valued function 
J(c[k]), called the error criterion , is a quantitative 
measure of the equalizer performance in reducing ISI. Such 
a function is also called an index of performance, cost 
function, and objective function. The ISI remains constant 
when the gradient VJ(c[k])=0, and from (2.2), no change 
is made in the tap-gains, i.e., the algorithm converges. 

Since the purpose of the algorithm is to determine a tap- 
gain vector c which minimizes the ISI, it is necessary to 

A 

define J(c) so that the algorithm converges to d. One way 
of doing this is to select J(c) to be convex. It then 
follows that any value c satisfying VJ(c) = 0_ also 
minimizes J(c) [29.]. 

Subject to these considerations, a convenient choice 
for the error criterion, analogous to the continuous integral- 
squared error [30], is the quadratic mean-squared error 
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(MSE) criterion given by 

oo 

J(o) = £ e m 2 , (2.3) 

m=-» 

where e - y_ - is the error corresponding to the 
transmitted symbol d . Using the complex convolution 
theorem (A. 10), integrating around the unit circle of the 
z-plane, and noting that S(z) = Y(z) - D(z), it follows that 

J(c) = £ £(z -1 ) £ (z) z” 1 dz 

= c T A c - 2 c T b + d, (2.4) 

where 

A = jL^VCzjV^z -1 ) z -1 dz = jij^X(z)F(z)F T (z - 1 )X(z -1 ) z _1 dz 

b = *.®[>(z) X(z) D(z -1 ) + F(z" 1 )X(z" 1 )D(z)] z _1 dz (2.5) 

and 

d = jiy f D(z) D(z -1 )z -1 dz. 

By letting z = e^ w it is easy to show that both the M-vector 
b and the scaler d are real. The matrix A in (2.5) is 
shown to be hermitian and positive semidefinite in B.l. To 
insure that A is also real, hence real symmetric, it is 
assumed that the filter-sections are FIR filters so that the 
output of F^(z) corresponding to input X(z) has the form 
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K-l 

V m (z) = E v mk z_k - (2 - 6 

k=0 

This assumption will be shown valid in (2.26) for the 

particular equalizer considered. Now from (1.5) it follows 

that any arbitrary element a of A is real since 

mn — 


K-l 

a mn = 2Wjf V z > = E v mk v nk' < 2 ' 7 > 

k=0 

Since A is positive semidefinite, J(c) is a convex function. 
Thus the gradient algorithm is applicable and the gradient 
is given by 

VJ (c) = 2 A c - 2 b. (2.8) 

The unique minimum value of ISI thus occurs for any tap- 

. ^ /s 

gain vector c satisfying VJ (c) = 0, or equivalently when 

A c = b . (2.9) 

2.3 THE TRANSVERSAL EQUALIZER 

Since A is positive semidefinite it may be singular, 
in which base there will not be a unique solution c for 
(2.9) . However, it is possible to select the filter-sections 
so that A is strictly positive definite, as will be shown. 

In particular, if F m (z) = z~ m for m=0 ,1, . . .M-l, and all 
filter-sections share the same shift register, then a 
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transversal equalizer results. To distinguish this case 
the quantities F(z) , V(z) , A, and b in (2.5) will be 
identified by a tilde, e.g., A. It is shown in B.2 that 
A is positive definite, implying that it is nonsingular. 
Thus A~*^ exists and the optimal tap-gain vector is found to 
be £ = A b. The tap-gain error vector e is now defined 
by 




( 2 . 10 ) 


and is given at iteration k by e[k] » c[k] - c. By adding 

"'—i 

and subtracting the term b A b, (2.4) can be reformulated 
as 

J(c[k] ) = d - b T A -1 b + e T [k] A e[k] . (2.11) 

T *“ 

The term e Ik] A e[k] is a positive definite quadratic form, 
which is zero only when e[k] - 0_. In this case (2.10) 

A 

implies that c[k] = c, so the irreducable error equals 

J . = d - b T A -1 b. (2.12) 

mm — — — 

The correctable error after k iterations is 

J Ik] = §. T t k J * e[k]. (2.13) 

Substituting (2.8) into (2.2) and using (2.10) gives the 
algorithm as 


clk+l] = c[k] - y A elk]. 

Repeated substitution into (2.10) with decreasing index k 
results in the "telescoping" form 
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k-1 

e[k] * lTT (I “ Y h 
1-0 


e [0] 


- (I - Y A)* e [0] . 


(2.14) 


Further simplification results from the observation that, 
since A is real symmetric, all eigenvalues are real and 
there exists a unitary transformation Q such that 
Q A Q - A, where the diagonal elements of 
A = diag(X'Q, Aj.,...A^_^) are the eigenvalues of A and the 
columns of Q ere the corresponding normalized eigenvectors 
{g m ) of A. Upon making the substitution A = Q il Q (2.14) 
becomes 

e [k] = Q (I - Y A) k Q T etO] , 
which when substituted into (2.13) gives 

M-l 

J cor tk] = E ^ T t°]S m ) 2 • < 2 - 15 > 

m=0 

It has been observed by Chang [25] that the 

eigenvalue spread of A, i.e., the difference between the 

largest and smallest eigenvalues, is a crucial factor in 

the rate of convergence. If all eigenvalues are equal 

then J can be reduced to zero on the first iteration 
cor 

by setting y - 1 . Unfortunately, since A depends on 

the input sequence to the equalizer, there is no general 
way to control the eigenvalues without changing the 
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equalizer structure. A logical approach is to investigate 
a more suitable form for the filter-sections {F m (z)}, and 
this will be done in the next section. It is interesting 
to note that by selecting a more general algorithm the 
correctable error can be reduced to zero on the first 
iteration using the transversal structure. This case is 
detailed in B.3, 

2.4 THE GENERALIZED TRANSVERSAL EQUALIZER 

Perhaps the simplest extention of the transversal 
equalizer is that generalized equalizer structure in 
which each filter-section is itself a transversal filter. 
Such an equalizer will be called a generalized transversal 
equalizer , and a typical filter-section is shown in Fig. 

2.3. This equalizer can be implemented so that all filter- 
sections share a common shift register. For convenience it 
will be assumed that all filter-sections are of equal 
length N, which will be strictly less than M. The resulting 
filter-sections are described by 

N-l 

F (z) = f z n , m = 0 , 1 , . . J4 -1 . (2.16) 

m I-* mn 

n=0 

They can be considered as devices which form linear 
combinations of those N values of the channel output 
sequence contained in the equalizer shift register, or 
equivalently, as the following linear time-invariant 
transformation on this sequence. 
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(2.17) 

The MxN matrix in (2.17) is now augmented by adding M-N 

zeros to each row, and the resulting MxM matrix is 

denoted F . Notice that column m of F corresponds to the 

coefficients of F (z) . This corresponds to physically 

m 

extending the shift-register length from N to M and 
assigning tap-gains of zero to the new tap-gain multipliers, 
and is done strictly for mathematical convenience. Upon 
recalling that for the transversal case 

t 

^ -1-2 ~M+1 T 

F(,z) = (1, z , z ,..,z ) , the relation 

m 

F(z) = P F(z) (2.18) 

T f- 

is seen to hold. Noting that V(z) = F F(z) X(z) 

= F T ]/(z} it follows from (2.5) that A = F T A F and 
T ^ 

b « F b. Substituting these expressions into (2.4) 
gives the new MSE expression 

J(c) = c A F 1 A F <2 - 2 c F b + d. 


(2.19) 
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Taking the gradient of (2.19) gives 

V J (c) = 2F T AFc-2F T b^2 (A c - b) , (2.20) 

A A 

Thus any vector c satisfying A c = b minimizes the ISI 

and is consequently an optimal tap-gain vector. It is 

crucial to note here that such a vector exists even though 

A is singular because, as proved in B.4, b is always in 

the range space of A. The optimal tap-gain vector in 

A -I -X . 

this case is given by c = A b, where A xs the 
pseudoinverse of A as defined in A. 3. 

T -I 

By adding and subtracting the term b A b in 
(2.20) and using the facts from A. 3 that A A 1 b = b and 
A A" 1 A - A, the error criterion can be reformulated in 
terms of the tap-gain error vector 

e=c-c=c-A*b. (2.21) 


This gives 


J(c) - d - b T A -1 b + e T A e. 


( 2 . 22 ) 


T -n T 

The term d - b A b Is independent of the tap -gain -vector 

A 

c and represents the minimum error J(c) - while the 

quadratic term again represents the correctable error 
Jcor* symmetry of A follows from that of A upon 

noting that, for arbitrary m and n, a ik = a ki implies 
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I 
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M-l 

a mn YL 
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II 
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M-l 

M-l 

n 

M 
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Therefore, the matrix Q having the normalized eigenvectors 

{q m } of A as columns forms a unitary transformation that 

diagonalizes A according to Q A Q - A = diag (X^ , , . . . 

with {X > being the set of eigenvalues of A. By sub- 
m — 

stituting the gradient (2.20) into the algorithm and 
recognizing e[k] at iteration k from (2.21) it follows 
that 

c[k+l] * c[k] - y A e [k] . (2.23) 

A 

Subtracting c from both sides of (2.23) gives 
e[k+l] = e[k] - y A e[k] , 

which reduces to the following expression in terms of e [0] , 
e [k+l] = Q (I - Y A) k Q T e [ 0 ] . 

Finally, J cor M expressed in terms of this result is 

J cor [kl = ( 2. T e(0]) T (I - Y A) k A(i -y A) k (Q T e[0]) 

= E ( ^ I0] X m (1 ^ V 2k - 

m=0 


(2.24) 
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Since F, and hence A, is singular then one or more of 
the eigenvalues will be zero and the corresponding 
terms of (2.24) are then zero, too. This observation 
will prove useful in the next section. 

If an input sequence {d^} of length D is trans- 
mitted over the digital channel of (2.1) then the channel 
output is given in the transform doSttSfin by 


X(z) 


G{z) D ( z ) 


G-l D-l 


22 E 


'n 


d z _{m+n) 


Rearranging the summations and letting k = m+n results in 
D+G-2 

X (z) = £ 

k-0 
L-l 



k=0 


where L = D + G - 1 is the length of the channel output 

sequence. The channel output sequence {x. m > is given by , 

the discrete convolution of the channel input {d } and 

m 

channel unit-pulse response {g^}, as seen from (2.25). 

An analogous derivation for each filter-section of length 
N gives the filter-section output 


E g n d k-n 
n=0 


-k 


(2.25) 


V z > 


L+N-2 


k 


E 


E 

VI — (1 



K-l 


A 



(2.26) 


f mn x k-n 
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where the sequence {v^} has length K»L+N-1. 

At this point it is necessary to make some 
assumptions which will lead to the selection of the 
filter-section length N and the number of such sections M. 

A 

First, a finite upper limit L of the channel output 
sequence length L is assumed known. Since a FIR channel 
was assumed in Section 2.1 and the transmitted initializing 

A 

sequence {d^} is known, L certainly exists. Second, it 
is assumed that M is at least as large as the length of 
the filter-section output sequences K, i.e., 

M ;> N + L - 1 . (2.27) 

A 

For channels with L very large it may be necessary to 
truncate the channel unit-pulse response in order to 
obtain a reasonable value for M. This has only a slight 
effect on performance when the truncated terms are small. 
This last condition (2.27) insures that the output sequence 
length can always be extended from K to M by appending M-K 
zeros to the transmitted sequence ' {d} . The selection of 
numerical values for M and N will be considered in the 
next section. 

The convolution relation (2.26) can now be written 
in the following matrix form 
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1 

o 

> e 

1 


Xq 0 ... 0 

v ml 

m 


X 1 x 0 
• » 

• 

r 


! x i 

V 

m,K-I 

" 

x l-i : 

• 


? X L-1 X 0 

« 


• 

. 0 k x, 

• • -L 



' . 


or more concisely as 


“mO 


’■ml 


'm,N-l 


v = X f , m = 0 f 1,...M-1. (2.28) 

— m — — m 

The MxN matrix X has rank N since the NxN principle minor 
is lower triangular with nonzero diagonal elements . The 
filter-section output vector V is an M-vector, and the 
filter-section tap-gain vector has length N. With 
this vector notation an arbitrary element a^ of A can 
be written using (2.7) as 


a = , V \ , (2.29) 

mn m — n / 

which is the correlation of the output sequences from 

filter-sections m and n, i.e., A is a correlation matrix. 

It can now be seen from (2.24) that J is reduced to 

zero on the first iteration when all nonzero eigenvalues 

of A are equal, since the convergence factor need only 

be set at y = L 1 for ^ 0. Conditions for which 
1 m m 

A = 0 or A = y ^ are established in the next section, 
m m 
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2.5 CONDITIONS FOR ONE -STEP INITIALIZATION 

Since the elements of A depend on the filter- 
section output vectors ' their relation to the 

eigenvalues of A is most important. An ideal solution 

A 

would be to determine a set {f_} of ideal filter-section 

—m 

tap-gain vectors that would make the corresponding set 

A 

{ v } orthonormal for a given channel output . It then 
— m 

follows from (2.29) that A is the identity matrix with all 
eigenvalues equal to unity. Unfortunately, (2.28) represents 
M scalar equations in N unknowns with N < M which cannot in 
general be solved exactly for f^. A "best" approximate 
solution does exist, however, and is given by 


f = X" 1 v , (2.30) 

— m — 

where X _I is the pseudoinverse of X as defined in A. 3, and 

f is the actual filter-section tap-gain vector. The 
— m 

actual filter-section output vector v and the ideal 


vector y_ Which it approximates can be related by 


V = X f = 
— m m 


-I 

X X y. . 



(2.31) 


From A. 3 (X X~^)^ X X 1 = XX 1 , and substitution of this 
result along with (2.31) into (2.29) gives 


mn 




*m' * * 


-I 

XX v 
-in 


> 


(2.32) 
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Using the orthonormal sets {^ 0 , and 

£.1'***£n- 1^ and the 9 eneralized spectral repre- 
sentations for X and X 1 as defined in A. 3, it is easy to 
show that 


X X 


-I 


N-l 


E s-C 


( 2 . 33 ) 


Furthermore, the orthonormal set {ij^} is linearly 
independent and thus constitutes a basis for the real 

A 

vector space of dimension M to which and y^ belong. 

A 

Consequently, can be represented in this basis by 


v_ = 
-m 


M-l 


n=0 


mn 


%n t 


( 2 . 34 ) 


where h^ i^^is the coordinate in the direction 

A 

of Now let H, V and Y be the matrices with columns 

A 

coinciding with the ordered sets • an< ^ ' 

respectively. Then making H an orthogonal matrix also 

A 

makes V orthogonal, as is easily shown by observing that 
¥ is unitary or equivalently, that 


<v !«>- E 


M-l M-l 

E h mk h ni<— k' U > 

k=0 i=0 


M-l 

E 

k=0 


^mk ^nk ^ 
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for all m and n. Substituting (2.33) and (2.34) into (2.32) 
gives an expression for any element a^ of A in terms of 
the set {Ji^} , 

N-l 

a = Y* h , h , . (2.35) 

mn §—* nut nK 

k=0 

To make the set {h } orthonormal the matrix H must be 

— m — 

orthogonal. Before discussing the selection of such a 
matrix, however, the eigenvalues of A will be determined. 

By defining the MxM matrix V as having columns 
coinciding with and letting P - XX 1 , (2.31) and 

(2.34) can be combined into the matrix equation 

A 

V=PV = P¥H . 


Furthermore, observing from (2.29) that A can be written 
as A = V T V, it follows that 

A = H T P T P ¥ H . (2.36) 

It is shown in B.5 that P represents a projection operator 
so that p' 1 P = P^ = P , and in B.6 that P equals the 
block diagonal matrix diag(I^ | 0^_^) . Since H is orthogonal, 
(2.36) implies that A is similar to diag(I^ I 2^-N^ thus 

has the same eigenvalues, namely, N equal to unity and M-N 
equal to zero. The desired conditions for one-step 
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initialization have now been found, as can be seen from 
setting y = 1 in (2.24) and observing that each term of 
the summation is then zero. 

There are an infinite number of orthogonal matrices 
which could be used for H, however the class of Hadamard 
matrices, discussed in A. 2, is particularly well suited 
for generating H. Such matrices are binary in the sense 
that all elements are either +1 or -1, have orthogonal 
columns, and are known to exist for all orders which are 
multiples of four over the range of feasible values of M. 
Consequently, normalizing an MxM Hadamard matrix by 1//M 
produces an orthogonal matrix with all elements equal in 
magnitude. This is particularly convenient for digital 
computation since the normalization can often be done 
near the end of the computation rather than on each element 
of A initially. 

The filter-section tap-gains can be computed 

-I 

without actually computing the pseudo inverse X , since 
substituting (2.33) and (2.34) into (2.30) gives 



h mk kt ' 


m = 0 , 1, . . . 


M-l. 
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This expression involves only the first N rows of H and 
the eigenvalues {y m 2 } and eigenvectors {^1 of the NxN real 
symmetric matrix X T X as discussed in A. 3. A number of 
numerical methods exist for finding eigenvalues and 
eigenvectors which are particularly simple for real 
symmetric matrices. Some of these methods are discussed in 

C. 3. 


2.6 SIMULATION RESULTS FOR NOISE-FREE CHANNELS 

The equalizer initialization theory developed in this 
chapter was tested by using a digital computer simulation. 

A detailed description of the simulation, including a 
complete listing, is found in App. C. Since this simulation 
was written for a maximum of 32 filter-sections, i.e., 

M £ 32, there is a practical limit to the length of the 
channel unit-pulse response as indicated by (2.27). Under 
this restriction, three digital channel models were selected 
for their representative characteristics. In the absence 
of an actual data channel, two models were taken from the 
literature and one was determined in the laboratory. 

The first channel was a low-pass filter (LPF) having 
a normalized cutoff frequency of 1.0 Hz and a 40 dB/decade 
roll-off. It is described mathematically by the Fourier 
transform pair 
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G 1 (f> = (l+jf) _2 -W g x (t) = 4 tt 2 te~ 2vt , t > 0. (2.38) 

By sampling at 0.1 s intervals and truncating all sample 
values which are less than 0.02, the unit-pulse response 
{2.11, 2.25, 1.80, 1.28, 0.85, 0.55, 0.34, 0.21, 0.12, 

0.07, 0.04, 0.02} of length 12 is obtained. Although 
this channel is rather artificial for actual data communi- 
cations, it is useful to demonstrate equalizer performance 
with nonnegative samples which are relatively large with 
respect to unity. It also serves for comparison with an 
ideal bandlimited channel and the previous work of Burlage 
and Houts [31], An experimentally convenient and widely 
used measure of ISI is the binary eye opening and the 
related peak-distortion criterion [32], Assuming that 
the second sample corresponds to the desired unity output, 
then the normalized peak distortion is 3.284, and thus 
considerably more than the value of 0,5 required for 
eye-closure. 

For the second channel model an electronic bandpass 
filter (BPF) [33] was used which had a maximally- flat 
passband with break frequencies of 200 Hz and 3.0 kHz 
and 60 dB/decade roll-off. The unit-pulse response of 
this channel was determined in the laboratory by exciting 
the filter with white noise and crosscorrelating the filter 
output with delayed versions of the white noise. Assuming 
a sampling interval of 50 ys , the unit-pulse of length 11 
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was {0.50, 1.15, 0.50, -0.25, -0.50, -0.35, -0.20, -0.10, 
-0.15, -0.10, -0.05}. This channel passband approximates 
the amplitude characteristic of a voice-grade telephone 
line, with zero phase at midband and tending to ± 270° 
in the stopbands. Since the normalized peak, distortion is 
2.348, the binary eye is again well closed. 

The third channel model was taken from a paper by 
Wescott [343 , and represents the unit-pulse response of a 
vestigial-sideband (VSB) amplitude-modulated data link 
composed of two carrier systems and 100 miles of loaded 
cable. When sampled at 0.4167 ms, corresponding to a 
transmission rate of 2400 baud, the unit-pulse response of 
length 9 was reported to be {-0.05, 0.05, -0.20, -0.05, 

0.90, 0.12, 0.15, 0.05, 0.03}. The peak distortion value of 

0. 778 indicates a closed eye once again. Unlike the previous 
channels, this one has a relatively large number of pre- 
cursors as well as a phase characteristic that is quite 
nonlinear . 

For an approximation to a noise-free channel the 

SNR was programmed for 100 dB as computed on a per-bit basis, 

2 

1. e., E^/N^ = l/o • A single unit pulse was transmitted 
and the MSE for this pulse was computed in accordance with 
the error criterion (2.3). Simulation runs were made to 
study the effects of the following parameters on equalizer 
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performance; the delay CD) between the arrival at the equa- 
lizer of the first precursor of the channel unit-pulse response 
and the occurrence of the locally generated reference pulse, 
the number of shift-register stages (N) in each filter-section, 
and the number of filter-sections (M) each of- which has an 
associated tap-gain multiplier c m » It will also prove con- 
venient to introduce the notation (M,N,D) to describe the 

equalizer defined by these three parameters. For the study 
of D two different sets" of values were used for the para- 
meters M and N. In order to obtain as wide a range of D as 
possible the relatively large values of 32 and 20 were 
initially used for M and N, respectively. The second set 

A 

of values was selected so that N » L,.the length of the 
unit-pulse response of the particular channel considered, 
and M was then set to the smallest multiple of four con- 

A 

sistent with the constraint (2.27) i.e.,M^N+L-l. 

Thus a comparison is obtained between relatively large and 
small equalizer structures. For the study of N, M was again 
set to 32 and D assigned a nominal value suitable for the 
channel considered. Selection of this value will be dis- 
cussed shortly. In studying M, N was assigned the largest 
value consistent with (2.27) and D was again set to some 
nominal value. Finally, the equalized unit-pulse response 
of all three channels was determined using (32,20,D) and.* 

A 

(# y L,D) equalizers with M the smallest feasible value and with 
D set to the corresponding optimal value. 



43 


Selection of the nominal value of D involves several 
considerations. The basic convolution relation describing 
the channel is obtained from (2.26), which for the noise- 
free case and filter-section m reduces to 

k k 

v mk “ £ f mn g k-n ” ^ f m,k-n g n' (2.39) 

n=0 n=0 

where f^ is defined to be zero if k > N, and {g^} is the 
channel unit-pulse response. If there are i precursors 
{g Q , then (2.39) implies that the main pulse 

corresponding to a binary one cannot affect the 
equalizer output before i shifts of the N-stage register. 

It is intuitive to expect the best equalization to occur 
when both precursors and postcursors are able to influence 
the equalizer output sample corresponding to the transmitted 
unit pulse. However, since the MSE (2.3) depends on the 
filter-section outputs at each of the N+L-l sample times 
required for the channel output sequence to pass completely 
through the equalizer, this may not always be true. In 
general the optimal delay can be expected to be somewhat 
larger than i, provided the shift-register length is on 
the order of 2i or greater. For shorter lengths the 
optimal delay may be less than i in order to retain a 
sufficient number of precursors in the shift register for 
the best equalizer output. Since the equalizer tap-gain 



settings depend on D and the output depends simultaneously 
on all tap-gain multiplier outputs, the filter-sections 
cannot be treated independently in determining the optimal 
D. Consequently, the optimal D must be determined by trial 
and error, which is not feasible to incorporate in the 
initializing algorithm considered here; therefore, a nominal 
delay must be selected for a given type of channel. 

Additional insights into the selection of the 

nominal delay can be gained by considering the frequency 

domain. Since each filter-section is a transversal digital 

filter with, in general, no symmetry in the set of multipliers 

{f }, the phase vs. frequency characteristic is nonlinear 
mn 

and the corresponding delay cannot be expected to be 
constant. Rather, different frequencies will be delayed 
by different amounts. This same nonlinear phase charac- 
teristic is also typical of analog communication channels 
where it produces delay distortion . For example, the 
voice-grade telephone lines commonly used in data trans- 
mission typically employ such phase-distorting devices as 
loading coils and bandpass filters. Delay distortion tends 
to cause relatively long precursors in addition to an 
average time delay of the signal maximum, and also may 
cause a skewing of the signal peak which results in either 
a small delay or advance of the signal maximum from the 
average time delay. 



45 


The transmitted signal also undergoes amplitude 
distortion in the channel. This effect is due to nonuniform 
amplitude vs. frequency characteristics in both the channel 
passband and stopbands, and it causes dispersion of the 
transmitted pulse with the resulting precursors and post- 
cursors. Based on these considerations the nominal value 
of D for the LPF channel was taken to be approximately 
0.25 N since the channel has a nearly linear phase 
characteristic over most of the passband. A nominal value 
of approximately 0.50 N was chosen for the VSB channel in 
view of the loading and modulation, whereas, for the BPF 
channel, where the degree of phase nonlinearity is somewhere 
between the other channels, simulation runs were made using 
each of the previous choices. 

Results of the study of parameters D, M and N for the 
three aforementioned channels are presented in Figs. 2.4, 2.5 
and 2.6. Corresponding numerical data is found in Tables 2.1 
and 2.2. The effect on MSE when the delay D is varied while 
holding M and N fixed is shown in Part (a) of each figure. 

The anticipated variation of delay with the channel phase 
is apparent, and it can also be seen that the nominal values 
of D are reasonable. 

It is interesting to note that for the LPF channel 
with M set to 32 any value of delay less than six 
essentially minimizes the MSE. Because of the large 
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amplitudes of the first two samples of the channel output, 
either could be considered to correspond to the transmitted 
unit pulse. Hence the channel delay is small and could 
easily be zero, depending on which sample the equalizer 
chooses to correspond to the unity output. Due to the 
nature of the particular filter-section tap-gains computed 
for this channel, the smallest residual MSE values occur 
when the channel output sequence begins to enter the shift 
register, rather than after most of the postcursors are 
included as might intuitively be expected. 

The variation in performance with shift-register 
length (N) is shown in Part (b) of Figs. 2.4 through 2.6. 
Larger values of N produce better performance, as can be 
clearly seen. For the BPF channel there is little difference 
in the results for nominal delays of approximately 0.25 N 
and 0.50 N so long as N is larger than the channel unit- 

A 

pulse response length CL) of 11. However, when N becomes 

A 

less than L, the value of 0.50 N is clearly superior. 

This is as predicted, since the short shift-register lengths 
cannot hold complete information on both precursors and 
postcursors, and a delay of 0.25 N results in a relatively 
large absence of postcursors. The improvement in performance 
with an increasing number of filter-sections CM) is also 
clearly seen from Part Cc) of each figure. 
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The unit-pulse response of the respective channel is 
shown in Part (d) of each figure both before equalization 
and for the smaller equalizer case. Thus the plotted unit- 
pulse response represents a relatively small equalizer 
structure. Although the resulting distortion is quite 
small in all cases, it can be further reduced by increasing 
the equalizer size. This can clearly be seen from Table 
2.3 which gives the equalized responses for both large 
and small size equalizers and the optimal values of D 
pertinent to each case. From the table it is also apparent 
that the largest deviation of the equalized response from 
the desired values is about 0.08 for the BPF channel, while 
the corresponding figures for the LPF and VSB channels are 
on the order of 0.005 and 0.015 respectively. 

The equalizer design developed in this chapter affords 
considerable reduction of ISI for the noise-free channels 
considered. All three channels were successfully equalized, 

A 

even when N was less than L,the channel unit-pulse response. 

This is an important observation because, as discussed in 

C.3, the amount of computation required is proportional to 
3 

N . The effects of background noise and channel parameter 
variations will be analyzed in Chapter 3. 
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Fig. 2.6 Noise-Free Equalizer Performance for VSB Channel 
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TABLE 2.1 

EFFECT OF DELAY (D) ON MSE FOR NOISE-FREE CHANNELS 
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TABLE 2.2 

EFFECT OF SHIFT-REGISTER LENGTH (N) 
AND NUMBER OF FILTER-SECTIONS (M) ON 
MSE FOR NOISE-FREE CHANNELS 


LPF 

CHANNEL 

BPF 

BPF 

VSB 

32 

32 

32 

32 

0.25N 

0.25N 

0.50N 

0.50N 


MSE 


(xlCf 6 ) 

--- 

<*io 3 ) - 

— 

>10^ 

322.90 

139.50 


67.23 

229.29 

108.80 


67.03 

124.50 

83.76 

49.84 

67.14 

107.30 

65.49 

48.16 

66.88 

57.50 

60.02 

11.19 

50.40 

52.71 

34.12 

9.25 

48.35 

26.77 

32.35 

1.54 

39.23 

22.37 

23.98 

1.28 

35.31 

16.83 

21.82 

0.30 


13.68 

20.31 

0.20 




0.05 




CHANNEL 




LPF 

BPF 

BPF 

VSB 

N 

M-ll 

M-10 

M-10 

M-8 

D 

0.25N 

0. 25N 

0.50N 

0. SON 

M 

(*10 -6 ) 

MSE 

-<*io“ 3 ) 



16 

67.07 

229.20 ' 

108.80 

49.85 

20 

1 50.10 

107.3 

65.52 

11.19 

24 

35.06 

52.71 

34.14 

1.52 

28 

24.85 

22.37 

23.97 

0.30 

32 

13.42 

13.68 

17.02 

0.05 












TABLE 2.3 
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UNIT-PULSE RESPONSE WITHtd^) AND WITHOUT (x^) EQUALIZATION 













CHAPTER 3 


ON-LINE OPERATION OF THE INITIALIZED 
EQUALIZER WITH NOISY CHANNELS 

In this chapter the generalized transversal 
equalizer design of Chapter 2 is modified so that near- 
optimal operation is maintained on noisy channels during 
the message transmission which occurs after initialization. 
Since it is impractical to compute the exact gradient 
VJ(c) of the error criteron (2.3) for noisy channels with 
ISI , a computationally simple method for estimating this 
quantity is presented. Some convergence properties of the 
resulting algorithm are then developed, and an eigenvalue 
condition is obtained which insures convergence in a 
stochastic sense under a stated set of assumptions. When 
these assumptions are violated the possibility of divergence 
exists, and a simple test is then proposed to determine 
when reinitialization is desirable. Finally, the effect of 
noise on both the initialization procedure and the on-line 
tracking is studied by means of a computer simulation using 
the same channel models as in Chapter 2. 
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3.1 THE NOISY ON-LINE ERROR CRITERION 

Under the noise-free hypothesis of Chapter 2 the 
equalizer initialization procedure used a single isolated 
initializing pulse and there was only one iteration 
of the algorithm (2.2). It was feasible to compute the 
gradient of the error criterion exactly, at least to within 
the limitations of numerical accuracy, since all nonzero 
error samples resulting from the initializing pulse 
were available separately for the summation resulting from 
(2.3), i . e . , 

oo 

ivj(c I o] ) =z: e k ve k • 

" k=-» 

However, when the equalizer is operated in a decision- 
directed mode with a noisy channel, then the situation is 
less straightforward. First, because of the channel noise 
the equalizer output sample is a discrete random variable. 
Consequently, the equalizer output and the error criterion 
are random as well, and some type of statistical averaging 
is indicated. Second, as the message pulses are no longer 
isolated at the channel output, but rather overlap to 
produce ISI, it is necessary to notationally distinguish 
between error samples corresponding to different transmitted 
pulses. Since the on-line equalizer will be adjusted once 
for each transmitted pulse, it will be convenient to 
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denote the error sample at iteration k by e[k], Notice 

that in general e [k] depends on the precursors and post- 

til 

cursors of neighboring pulses as well as the k pulse. 

2 

By taking the expected value of c Ik} the error criterion 
X{c[k]) for the noisy on-line case, corresponding to (2.3), 
is given by 

J(c[k]) = E{e 2 [k]}, (3.1) 

Where E is the mathematical expectation operator. 

The on-line algorithm for adjusting the equalizer 
tap-gains now becomes 

c[k+l] = o[k] - 1 y VJ (c (k] ) , (3.2) 

where y is again the scalar convergence factor, and the 
gradient expression is 

| VJ(c[k]) = VE{ e 2 [k]}. (3.3) 

To simplify (3.3) the filter-section output vector v[k] is 

defined to be v[k] = (v^[k], v 1 [k] , . . . v M _^ [k] ) T where the 

component v^Ik] is the output value of the m filter- 

section at iteration k. Notice that v[k] is distinctly 

different from v defined in (2.28). If the decision- 
— m .. 

directed estimate d[k] of the k transmitted message pulse 

A 

is denoted by d[k] , then the error e[k] is given by 
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e[k] = y{k] - d[k] = <^£ [k] , v[k£>- d[k], 
and upon substitution into (3.3) there results 

ivE{e 2 tkJ> = E{ e [k] v [k] } 

= E{v[k) v T [k] c [k] } - E{d [k] v [k] } . (3.4) 

Although both variables e [k] and v [k] are explicitly 
available in the equalizer and their product is easily 
calculated, taking the expectation of the result on-line 
is not practical. To circumvent this difficulty the 
gradient (3.4) will be estimated using the so-called 
"noisy" estimate method suggested by Widrow [5]. The 
resulting expression is 

I VJ (o Ik] ) = e [k] v[k], (3.5) 

which is an unbiased estimate since 

E{ e Ik] v[k] } =ivj(c[k]). 

The algorithm (3.2) is now approximated by the easily 
implementable form 

cCk+1] - c[k] - y e [k] v[k] . (3.6) 

Equation (3.6) involves no explicit a priori 
statistical information, but it is based on the decision- 
directed hypothesis that the binary sequence of equalizer 
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output decisions closely approximates the transmitted 
message sequence. It should also be noted that the 
initialization procedure requires a priori information in 
the form of the locally-generated initializing pulse replica. 
When a priori statistical information is available during 
on-line operation more accurate algorithms can be obtained, 
although the gradient determination generally involves more 
computation. Several such cases are discussed in B.7. 

The equalizer structure developed in Chapter 2 can 
now be modified for on-line operation using (3.6). Fig. 3^1 
shows the complete equalizer block diagram. All variables 
shown correspond to the output sample y[k], i.e. , the 

i 

delay caused by the channel and filter-sections is suppressed 
in the notation. , The channel noise n[k] is added to the 
deterministic channel output s [k] to obtain the equalizer 
input x[k]. A filter-section output v m [k] is the weighted 
sum of those N input samples located in the shift register 
when y[k] occurs, and the weights are the filter-section 
tap-gains given by (2.37). The equalizer output is now 
given in vector form by 

y [k] = ^ c[k] , v[k]^> 

The shift to on-line operation involves only two 

/ 

switchings as shown in Fig. 3.1. First, the equalizer 
output decisions rather than the reference generator 




Fig. 3.1. Equalizer Structure With On-Line Tracking 
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output are fed back to produce the error signal in a 
decision-directed manner. Second, since the on-line 
computations are considerably simpler than those during 
initialization, it may be advantageous to use a set of 
digital multipliers to compute the on-line gradient 
estimate in lieu of the adaptive computer used for 
initialization. This is particularly true when the 
initialization computation is done by a time-shared general 
purpose computer. Of course, if the adaptive computer is 
continuously available the second switching operation is 
unnecessary. 

A final design consideration involves the detection 

of erroneous operation. One simple method for doing this 

is to compare | e [k] | with some predetermined positive 

threshold value e and to reinitialize the equalizer 

max 

whenever e is exceeded a given number of times in a 
max 

specified period. 

3.2 CONVERGENCE PROPERTIES OF THE ON-LINE ALGORITHM 

Before investigating the convergence properties of 
the on-line algorithm (3.6), several preliminary con- 
siderations are necessary. Since the noise. £s random, this 
convergence must be considered in some statistical sense. 

It is also mathematically expedient to model the noise as 
a stationary random process with zero mean and finite 
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variance. Although some forms of background noise are not 
stationary, e.g., switching and impulse noise, other forms 
such as thermal noise do fit this model reasonably well, 
particularly if the time frame is not too long. Furthermore, 
if the equalizer setting deviates too much from the optimal 
setting, then the excessive-error test on the resulting 
errors elk] will indicate the need for reinitialization. 

Taking the expectation of (3.6) under the stated 
hypothesis and assuming c[k] is statistically independent 
of v[k] gives 

c [k+1] = c [k] - Y (R vv c Ik] - rj v ), (3.7) 

where R = R vy (0) = E{v[k] v T Ik]} is the filter-section 

A ^ 

correlation matrix , r dy (n) = E{d[k] v[k] } defines a cross- 
correlation vector, and c[k] now denotes the expectation of 
the tap-gain vector. A channel output M-vector x[k] at 
iteration k consists of the current shift-register contents, 


which produce y[k], for the first N components, and the most 
recent M-N samples to have passed through the register as 


the remaining components. By now defining the crosscorrela- 
tion matrix R r = R v (n) = E{x[k]x T [k]} and vector r, = 

— XX "XX “ 

E{dlk]x[k]}, and noting that v[k]=F T x[k) , it follows that R vv = 



I and £dv 


T 

F r^. Under the mild assumption that the 


signal s [k] and noise n[k] are independent, it is shown in 


B.8 that R is positive definite and thus nonsingular. It 

“XX 

can then be shown by the same procedure as in B.4 that r^ v 
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lies in the range space of R , and consequently (3.7) is 

^ “I 

convergent for c = R w r^ v . Although this is not, of 
course, the exact optimal solution of the original 
deterministic algorithm (2.2), it is the optimal expected 
solution of the on-line algorithm (3.6) in which the 
"noisy" gradient estimates are used. 

To study the convergence of (3.7) c is first sub- 
tracted from both sides of the equation which is then 

A 

written in terms of the tap-gain error vector e[k] = c[k] - c 
as 

a [k] = (I - y R vv ) k_1 ell]. (3.8) 

Since R vv is a correlation matrix it is symmetric and 
positive semidef inite. Consequently, all eigenvalues {p^} 
are both real and nonnegative and there exists an ortho- 
normal set of eigenvectors Letting 0 and 0 denote the 

index sets corresponding to the positive and zero eigenvalues 
of R , respectively, a unique orthogonal decomposition 
e[k] * e + [k] + e°[k] exists, so that<^e + [k], e°[k]^> = 0. 

In this decomposition e[l] is written as 

e [1] = e + [l] + e°tl] = £ + £ “jli (3 * 9) 

ie© ieO 

where {a^} are the appropriate coefficients for the basis 
{£^}, Substitution of (3.9) into (3.8) now gives 
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e[k] = e + [k] + e°[k] = 


E 

ie0 


a i (l-YP i ) 1C_1 £j_ + 


E a i— i . 


ie0 


(3.10) 


The last term of (3.10) is simply e°Il] as can be seen from 
(3.9), i.e., e°[k) will converge if the convergence factor 
y is chosen to insure that 


iee { U - y Pi I > < (3- 11 ) 

since in this case 

l|e + Ik]||<^ (H-TpJJIIE “i^ i ll < ll- + • 

i£0 

Therefore, the condition to insure convergence of (3.7) 
is found from (3.11) to be 


Y 


< 2/p 


max 


r 


where 


p > p. for all i. 

'max l 

3.3 SIMULATION RESULTS FOR NOISY CHANNELS 

To study the combined theory for equalizer initiali- 
zation and on-line operation, the computer simulation used 
in Chapter 2 was expanded to reflect the theoretical results 
of this chapter. In particular, provisions were made for 
generating a random binary sequence to simulate the trans- 
mitted message, for adding gaussian noise with specified 
statistics, for adjusting the equalizer tap-gains at each 
bit time, and for detecting erroneous operation. The 
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computer simulation is described in some detail in C.l and 
C . 2 . 

The equalizer performance was again evaluated over 
the three channels discussed in 2.6, i.e. , the LPF , BPF , 
and VSB channels, but the SNR was varied from the essen- 
tially noise-free value of 100 dB down to 15 dB. At 15 dB 
the noise standard deviation is 0.1778 and the probability 
of the noise alone causing a bit error is 0.0025. When 
this effect is superimposed on the distortion due to ISI, 
it is apparent that the equalizer alone should not be 
expected to perform well when the SNR is in this range. 

In practice such low SNR situations are commonly handled 
by coding, and are usually incompatible with the high 
baud rate to bandwidth ratios that justify sophisticated 
equalization. 

In addition to the equalizer parameter set (M,N,D) 
defined in 2.6, the convergence factor y must be specified 
for the on-line algorithm. It was found by trial and error 
that a good selection was 0.Q2, which represents 4% of the 
decision threshold value. Larger values of y tend to 
increase the bit error rate when there is noise in the 
channel, and can lead to divergence if the increase is too 
large. On the other hand, smaller values result in slower 
equalizer convergence rates. The transmitted message 
was simulated by repeatedly generating a 255 bit pseudo- 
noise (PN) sequence, and the error criterion was taken to 
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2 

be the average MSE e [k] after 500 iterations of the on-line 
algorithm. In determining e [k] the transmitted bit d[k] was 

A 

used in place of the output device decision dlk] for pro- 
gramming convenience. This substitution has no effect when 
there are no bit errors, and the effect is considered 
negligible for the SNR range considered. Data in which 
one or more isolated bit errors occurred is indicated by 
an asterisk in the appropriate figures. 

The first parameter studied was the delay (D) 
between the arrival at the equalizer of the first precursor 
of the channel unit-pulse response and the occurrence of 
the locally-generated reference pulse for initialization. 
Since D depends on such nonlinear effects as the channel 
and filter-section phase characteristics, it is impractical 
to determine the optimal value analytically. As in 2.6, 
simulation runs were made for a relatively large equalizer 
structure, (32,2Q,D), and for one for which the equalizer 
shift-register length N equaled the channel unit-pulse 
response, i.e., the structure is (M,£,D). The simulation 
results are shown in Part (a) of Figs. 3.2, 3.3, and 3.4, 
from which the effect of D as well as channel SNR at 20 dB 
and 40 dB can be seen. All three channels were equalized 
successfully for a SNR of 30 dB or greater using the larger 
equalizer structure. However, the smaller structure was 
unable to equalize the BPF channel, and the larger structure 
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was ineffective on this channel at 20 dB; consequently, a 
value of 30 dB vice 20 dB was used. The curves in Figs. 3.2 
and 3.4 indicate that there is no significant improvement 
in performance for the larger structures provided D is 
properly chosen, and in fact for the VSB channel the smaller 
structure performs slightly better. However, there is a 
much wider range of feasible values of D for the larger 
structures. Thus the possibility of poor equalization due 
to an inaccurate choice of D with a small structure must 
be weighed against the complexity and cost of realizing a 
larger structure. Since the accurate selection of a 
nominal value of D depends on knowledge of the channel, a 
larger structure is particularly suitable when the channel 
is relatively unknown, or may be expected to change 
significantly after initialization. Based on the con- 
sideration of 2.6 and the experimental results, nominal 
values for D of 0.50 N, 0.25 N, and 0.80 N were selected 
for the LPF, BPF, and VSB channels, respectively. 

The next parameter studied was the shift-register 
length (N) . The result of varying N with M set to 32 and 
D to the corresponding nominal value is shown for each 
channel in Part (b) of Figs. 3.2, 3.3, and 3.4. For high 
SNK values performance tends to either improve or remain 
nearly constant with increasing N. However, with the LPF 
channel the performance degrades slightly with N > 8. 
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Fig. 3.2 Equalizer Performance for Noisy LPF Channel 
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This may be attributed in part to the increased number of 
arithmetic operations, each of which is reduced in accuracy 
by the noise, that is inherent in the larger structures. 

It is also evident for the LPF and VSB channels that, for 
SNR £ 30 dB, there is a lower value of N above which there 
is little change in performance. As previously noted, 
this is important since the amount of computation required 
for initialization is reduced considerably with N. Again, 
however, it should be kept in mind that small values of N 
make the selection of D more crucial. It is interesting to 
note that for the LPF channel equalization was achieved 
when the SNR was 20 dB with a shift-register having only 
three taps. 

Part (c) of Figs. 3.2, 3.3, and 3.4 shows the effect 
of SNR on equalizer performance when both a single 
initializing pulse and an average of four such pulses were 
used by the equalizer to identify the channel during 
initialization. It is clear in each case that the use of 
multiple initializing pulses is quite beneficial where the 
SNR is low and the noise represents a significant percen- 
tage of the unit-pulse amplitude. This is not the case at 
SNR above 30 dB, so little advantage results from using 
multiple initializing pulses. 

i 

The LPF and VSB channels were successfully equalized 
with a 20 dB SNR? however, 30 dB was needed before the BPF 
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channel was equalized without bit errors. As can be seen 
from a comparison of the channel unit-pulse responses shown 
in Part (d) of Figs. 3.2 , 3.3, and 3.4, the BPF channel 
has more overshoot than the other two, and also has a 
number of large, negative postcursors. Thus it is expected 
to be the more difficult channel to equalize. The fact 
that no value of N was reached above which the equalizer 
performance remained fairly constant suggests that better 
equalization is attainable for larger equalizer parameters 
than the maximum of (32,22,D) permitted by the particular 
computer simulation used. There are a number of tradeoffs 
in selecting the equalizer parameters which can now be 
summarized. First, selecting large values of M and N gives 
a larger range of values from which to select the nominal 
value of D, which is particularly important if little is 
known about the channel a priori. On the other hand, 
smaller values of N greatly reduce the amount and time of 
computation required for initialization, while smaller 
values of M reduce the number of on-line multiplications 
required at each algorithm iteration. These considerations 
can also be thought of as a tradeoff between the accuracy 
desired from the equalizer and the cost that must be paid 
to achieve it. 



CHAPTER 4 


CONCLUSIONS AND RECOMMENDATIONS 

The design of a linear digital equalizer for reducing 
ISI on digital communication channels has been developed 
and evaluated. This design includes an initialization 
algorithm employing a digital computer, as described in 
Chapter 2, and a simplified on-line algorithm, developed 
in Chapter 3, which requires only a set of digital multi- 
pliers. The effects of the main equalizer parameters 
(M,N,D) and of channel SNR were studied by means of a 
digital computer simulation for three baseband channels. 

All channels were equalized successfully for a SNR of 
30 dB or greater, and two were equalized at 20 dB. 
Considering a digital computer as a sophisticated digital 
filter, the equalizer simulation is quite realistic except 
for the relatively long word length of 36 bits. 

There are a number of considerations to be taken 
into account when realizing an equalizer of the type 
presented here. First, the maximum shift-register length 
N is clearly limited by the digital computer used for the 
initialization computation. The amount and speed of 
digital hardware available for the shift register, 
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multipliers, and summers limits N, and M as well. Another 
important consideration is the determination of a nominal 
value for the delay D. Good estimates of D will generally 
require some a priori knowledge of the channel character- 
istics or some preliminary experimental data. An alter- 
native procedure, which might be necessary with a completely 
unknown channel, is to iterate the initialization procedure 
several times using different values for D. If the nominal 
value of D is uncertain, or if considerable variation is 
expected in the channel unit-pulse response, then a large 
value of N will make the equalizer performance less 
sensitive to D. The nominal value of D ranged from 0.25N 
to 0.8N depending on the choice of channel model. A third 
design consideration is the test for excessive error 
indicated in Fig. 3.1. Since single bit errors are highly 
probable when the SNR is low or the channel badly equalized, 
some type of averaging is desirable to avoid an excessive 
number of reinitialization commands. This could be done, 
for example, by keeping a running count of the number of 
times I e Ik] I exceeds a predetermined limit e over the 
last m iterations and reinitializing only if it is 
exceeded n < m times , where m and n are known integers . 

Such a scheme would also prevent short bursts of impulse 
or switching noise from causing a needless reinitialization. 
This scheme could be easily implemented with digital 
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hardware, as could more sophisticated schemes which might 
be devised for particular installations. 

There are a number of interesting areas which could 
provide motivation for further investigation. One 
particularly germane area is the equalizer performance on 
time-varying channels and the effect of the convergence 
factor y on this tracking. The computer simulation of 
App. C is easily modified for such a study. A second more 
theoretical area involves reduction of the amount of 
computation required for initialization. More efficient 
algorithms for finding eigenvalues and eigenvectors could 
perhaps be devised for particular equalizer installations. 
The use of a suboptimal procedure to estimate the initial 
equalizer settings might also reduce the initialization 
time considerably . The third area to be suggested concerns 
the effects of short word lengths and the related problems 
of truncation and round-off error. There has been some 
interesting theoretical work in this area 113] , and 
experimental simulations could be done using a minicomputer 
for which a typical word length is 16 bits. A fourth area 
for further investigation is the possibility of further 
simplifying the on-line algorithm by using only polarity 
information of some of the variables. This type of scheme 
has been called ’'hybrid" by Hirsch and Wolf [10] who studied 
it for transversal equalizers. 
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A final comment concerns the feasibility of physically 
realizing an equalizer of the type developed herein. Al- 
though the equalizer requires a considerable amount of 
computing equipment for ini ti aliz ation , the constant 
improvement in integrated circuit technology and the 
development of faster, more complex components make the 
production of such equipment appear feasible. Furthermore, 
the proliferation of time-sharing systems and inexpensive 
minicomputers makes large amounts of computing facility 
widely available. Thus the equalizer design presented 
herein is both theoretically and physically feasible. 



APPENDIX A 


MATHEMATICAL FUNDAMENTALS 

This appendix will provide both review and 
supplemental mathematical material. It deals first with 
selected topics from linear algebra and operator theory 
on finite dimensional inner product spaces. This theory 
is then specialized to matrices on real and complex vector 
spaces, with the subject of pseudoinverse matrices covered 
in a separate section. Finally, the classical z-transform 
is briefly discussed and some useful properties are 
derived. 

A . 1 LINEAR VECTOR SPACES AND OPERATORS 

A mathematical field F consists of a set of elements 
called scalars which satisfy a particular set of axioms [35]. 
The only fields used herein are the sets of all real and 
all complex numbers. A linear vector space over a field F 
is denoted (tt,F) or just U when F is clearly understood. 

The elements of U , ctalled vectors , also satisfy a certain 
set of axioms [ 36 ] and are related to each other and to 
the scalars by the operations of vector addition and scalar 
multiplication. Linear vector spaces are also called 
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vector spaces or linear spaces . Vector spaces are called 
real or complex according to whether the field is real or 
complex, respectively. A set of vectors {x^ X2 '****k^ i n 
U spans U if every vector £ in U can be expressed as a 
linear combination of the set, i.-e., if there exists a set 
of scalars {a^, such that 

K 

* - L \ . 

k=l 

A set {x^} is linearly independent if the only way that their 
linear combination can be made to equal zero is for all of 
the scalar multipliers to be zero. Any linearly independent 
set that spans U is called a basis of U, and the number of 
elements in any basis is called the dimension of U. Only 
vector spaces of finite dimension will be considered here. 

Let (t/, F) and (V, F) be vector spaces with a, 3 being 
arbitrary scalars in F and w, x, £ arbitrary vectors in U. 
Suppose that to each element of U there is assigned a unique 
element of V, then the collection of these assignments is 
called a map T from U into (/. Maps are also called functions , 
transformations , or operators , and these terms will be used 
interchangeably . 

The identity operator I is defined by I x = x, and a 
map T” 1 satisfying the relation T" 1 ! » TT ^ = 1 is called the 
inverse of T. The inverse operator may not exist; however, 
when it does then T is said to be invertible . The set of 
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scalars X for which the map (T - \ I) is not invertible is 
called the spectrum of T and will be denoted a[T]. When T 
satisfies the relation 

T (a x + 3 y) = a T (x) + 3 T (y) 
it is called a linear operator . It is interesting to note 
that a field is a vector space, and also the set of all 
linear operators from U into 1/ forms a vector space in which 
each linear operator is itself a vector ^pace [37] I 

There are many important maps between vector spaces. 
An i nner or scalar product is a scalar-valued function of 
two vectors, denoted ^ x, , which is defined by the 

following three axioms; 

1 • <(z' * y 

ii. x + 3 Y_r w > --< X, w y + s <£, w y 

iii. ^x, x ^ >0 with equality only when x = £ . 

A vector space with an inner product defined on it is 
called an inner product space . The inner product determines 
the relative directions of vectors, and two vectors are 
orthogonal if ^x, y ^ = 0. The notion of distance in 
vector spaces is indicated by a nonnegative scalar-valued 
map called a norm , |] x [| , which is also defined 

axiomatically by 

If « 2 II = M • || 5 |[ 

iJ -- IU + n (I 6 tun + nn n 

iii. || x || >0 with equality only when x = £. 
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A set of vectors {x_} is an orthonormal set if |(x„ |[ - 1 

— m 511 

for any in and any two distinct vectors are orthogonal . A 
vector space with a norm defined on it is a normed vector 
space . For the finite dimensional vector spaces considered 

herein it will be convenient to define the inner product of 

T . T 

two vectors, say x = (x^ x 2 ,...x K ) and y = (y^r Y2'** ,y K' » 

by 


K 


<*' £> A E x k *k* = 


* 

* Z , 


k=l 


which then defines a norm of x by 

II 2 . II 4 */<*' *> . 

A map 7, not necessarily linear, is called convex if 
for every value of a satisfying O^a-1, it has the following 
property; 

T( [1 - a] x + a y) ^ (1 - a) T (x) + a T (y) . 

If the inequality is strict, then T is strictly convex . 

For every linear operator T on an inner product 

A 

space U there exists another operator T called the 
adjoint of T defined by 

If T “ T A then T is called a self-adjoint operator , and if 
T T then T is a normal operator . A projection 

operator is one which is both self-adjoint and idempotent , 
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i.e., T T a T. Any map for which the adjoint equals the 
inverse is a unitary operator . 

A final important notion deals with the geometric 
interpretation of a map 7 from U to (/• The set of elements 
x in LI such that T x = £ is called the null space , or some- 
times the kernel, of T and is denoted N(T).' The set of 
elements y in t/ which can be written y = T x for some x in 
U is called the range space , or image, of T and is denoted 
R(T), The null space of T is related to the range space 
of its adjoint T A in the following way; every element 
in W(T) is orthogonal to every element of R(T A ). To show 
this note first that x belongs to W{T) if and only if 
<x > * > - 0 for any y in U, and since ^y, Tx^ = 

^T A x> y^then x is orthogonal to every w in U such that 
& 

w « T y for some y. But this set of w's is just the range 

a 

space of T , hence the desired result is proven. 

A useful consequence of this last result is that the 
A A 

range spaces R(T ) and R(T T) are identical. This can be 

a 

proved indirectly by first noting that if x is in N(T ) 
then (T T A ) £ ■ T (T^x) » T (0) = 0 so x is also in W(T A T) , 
Conversely, if x is in W(T A T) then |( Tx [| ^ = ^Tx, Tx ^ ® 
<*, T* r *>=<*,0> = 0, so that Tx = £ and x is in W(T). 

The proof is completed by noting that every element of R(T A ) 

a 

and R(T T) is orthogonal to every element of W(T) and 
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N(T ‘k T) , respectively* • Hence R(T^) and R(T^ T) must then be 
identical. 

A. 2 MATRICES 

An M X N matrix T is a rectangular array of elements 

having M rows and N columns? the element in row m and column 

n is denoted by t and the entire matrix T is sometimes - 

mn 

denoted by [t^] . By defining the sum of two matrices and 
the multiplication of a matrix by a scalar in the usual way , 
it can be shown that the set of all M x N matrices, with 
elements in the same field F, forms a vector space [ 33 ] . 

The way that matrices are used to represent linear 
operators is made precise by the following Repr e s ent at ion 
Theorem [393? 

Let (U, F) and (!/, F) be real or complex vector 
spaces with bases {u^» ^—1' — 

respectively, and let T be a linear map from U 
to i/. Then, with respect to these bases, T is 
represented by the M x N matrix T « [t^] where 
the elements of column k of T are the components 
of T(u k ) with respect to the basis {v^h 

When M = N then T is a square matrix and U is identical 
to i f, i.e., T maps U into itself. If M = N and 

u 2 , . - • u^ ) is another basis for U in which ¥ represents 
T, then there exists a nonsingular square matrix P that 
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relates these representations according to T = t P. 

Such a mapping is called a similarity transformation , and 

the matrices T and T are said to be similar . Thus similar 

matrices represent the same operator in different bases. 

For a square matrix T of order N the spectrum 

c. [T] is that set of scalars for which (T - X I) fails to 

be invertable, or equivalently, for which the determinant 

det(T - X 1) =0. This last equation is a polynomial in 

X of order N, so that a [T] contains exactly N elements 

which may not be distinct. These elements are called the 

eigenvalues of T, or sometimes proper or characteristic 

values. The spectrum is a property of the abstract 

operator T and does not depend on the particular repre- 

/*✓ 

sentation, e.g., T or T. This can be seen by recalling 
that T = P” 1 TP so that 

det (T- X I) = det (pf 1 T P - X P _1 P) 

= det (P -1 ) det (P) det (T- X I) = det(T- X I). 

Any nonzero vector. £ such that (T - X I) £ = 0, with X in 
a [T] , is called an eigenvector of T corresponding to the 
eigenvalue X. 

There are a number of special matrices that will 

prove useful herein, some of which represent the various 

operators defined in Section A.l. The square matrix 

a 

representing the adjoint T of an abstract linear trans- 
formation T on vector space li is the complex conjugate 
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H • 

transpose of T denoted T as can be seen from the following, 

\X’ !Z/=x CT x\* = X 1 T y - < T x, £ y . 

When T is self-adjoint T = T** and T is called a hermitian 
matrix . If the vector space U is real rather than complex, 
then T = T and T being self-adjoint implies T = T . Such 
a matrix is called real symmetric . The eigenvalues of an 
N x N hermitian matrix T are always real, and there always 
exists a set of N linearly independent eigenvectors 

which are orthogonal [4Q.].. In fact, since when £ is an 
eigenvector then £/ II £ II is too, then the set {p^} can 
always be taken orthonormal. Writing the eigenvector equator; 
Tp m = for m ^ 1 ' 2,...N in matrix form gives 

T * E-2 9 * * * En ^ = 9 &2 1 * * * Ejj ^ diag (X^, X 2 * • • • Xjg ) 

or, more concisely in terms of the matrices P and A, 

T P = PA. 

Furthermore, since ^p^, p^ ^ * 6^ and a similar relation 

— 1 H 

exists for the rows of P [41], then p_ = which implies 
P is a unitary matrix . Consequently, the similarity 
transformation A = T P that diagonalizes T becomes 

A = P H T P f (A. 1 ) 

which is called a unitary transformation . Hermitian 
matrices can be easily decomposed into their spectral 
representations by using unitary transformations in the 
following way. 
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N 

1= E a p h = £ • (A - 2) 

m=l 

When the vector space is real and thus T is real symmetric, 
the complex conjugate transpose is simply the transpose. 
Results analogous to the hermitian case exist for the 
general case where T does not have N linearly independent 
eigenvectors, but the theory is considerably more involved 
[42] . 

A projection operator is represented by a projection 
matrix T which is hermitian and has the idempotent property 
T 2 = T. It can be shown [43] that all eigenvalues are 
either zero or unity, and that the eigenvectors of the 
unity eigenvalues span the subspace into which the operator 
projects. Thus anyv projection operator is similar to the 
block diagonal matrix diag(lJO) = diag (1 ,1 , . . .1 , 0,...0) 
where the square matrices I and 0_ are of orders equal to 
the number of unity and zero eigenvalues, respectively. 

A quadratic form is a real-valued function of N real 
variables {x^} defined by 

N 

F(Xi, x 2 / • • • ^) = 

m=l 


N 


'mn 


*m 


n 


(A. 3) 


n=l 


which can also be written as 
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N N 

E r mn x m 2 + E E tr mn + r nm ) x m V 

m=l m=l n<m 

The value of this last expression remains unchanged if the 
substitution = 1/2' (r^ + r nm ) is made. By defining 
the real symmetric matrix Q 419^3 and letting x = (x^, 
x 2 ,...x n ) T , (A. 3) can be written in matrix form as 

F (x) = x T Q x = ^x, Q x ^ . 

In the more general complex case when Q is hermitian the 
function is still a quadratic form since then 

N N 

F <*> = E x m 2 + E E (<J mn + O x m x n / 

m=l m=l n<m 

* 

which is real due to q m + = 2 Refq^). 

A real symmetric matrix Q is positive semidefinite 
if, for any choice of x, the quadratic form F(x) = ^x, Q x^ 
is nonnegative. If F(x) = 0 only when x = 0_ then Q is 
positive definite . The eigenvalues of Q are related to 
the definiteness in the following way? if Q is positive 
semidefinite all eigenvalues are nonnegative, while all are 
strictly positive when Q is positive definite. If Q is 
positive definite it is invertible. 

The final type of special matrix to be considered 
in this section is the Hadamard matrix > which has orthogonal 


rows and columns, and is binary in the sense that all 
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elements are either 1 or -1. Except for the 2x2 case 



Hadamard matrices may exist only for orders which are 
integer multiples of four, and have been constructed for 
all such orders up to 100 [443. For those orders which are 
also powers of two, the construction is done recursively by 
starting with and using the iterative relation 

— k+1 A 

The construction of other orders involves the determination 
of quadratic residues [%5] and is thus somewhat more involved 
in theory but easily implemented using a digital computer 
as discussed in Section C.2. 



A. 3 THE PSEUDO INVERSE MATRIX 

Before defining the pseudoinverse matrix a number of 

preliminary notions are necessary. Let the M x N matrix X 

represent a linear map from a vector space (I to a normed 

vector space V of dimensions N and M, respectively. Assume 

that the field is the set of real numbers and that M £ N. 

T 

Then all eigenvalues of the real symmetric matrices X x 
T 

and XX ard £eal and nonnegative. The real property follows 
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T 

from symmetry; to show that X X has nonnegative eigenvalues, 
let v = X u for some arbitrary vector u in (I. Then 

0 < [ | v 1 [ 2 = v T v = u T X T X u f 

so X X is positive semidefinite from which the desired 

T 

result , follows . An analogous procedure holds for XX . If 

the rank of X is R < N, then it can be shown [46J 

that both X T X and XX T have exactly R identical nonzero 

eigenvalues, say {y 1 / y 2 /•••Pr K with all other 

eigenvalues zero. Furthermore, from symmetry it follows 
T T 

that X X and XX each have a set of orthonormal eigenvectors 

which will be denoted by the sets and 

&2 • • * • lfe.24^ ' respectively. These sets are ordered so 

2 

that and ^ correspond to y^ for m « 1, 2,...R. The 
rectangular matrix X can now be expressed [47] in the 
generalized spectral representation 

R 

* = E ik ik T • (A - 4) 

k-1 

When X is real symmetric, hence square^then (A. 4) reduces 

to (A. 2). This cart be shown by noting that XX T = X T X « X 2 

and, from the Frobenius Theorem the nonzero eigenvalues of 

X are {y 1 , y 2 ,...y R }. If £ is a normalized eigenvector of 

2 

X corresponding to y, it is also an eigenvector of X 

2 

corresponding to y . To show this note that (x - y I) E = 2. 
implies that (X 2 - y 2 X)£ * X (X £) - y X £ = X [ (X - y I ) £] 
= 0. Consequently, for this case {g^}, {^}, and {3^} are 
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all equal and the desired result follows. 

The N x N pseudo in verse matrix x ” 1 of X is now defined 
by 



R 


E 



(A. 5) 


One important application of this definition can be seen 
by considering the matrix equation X u = v. If v belongs 

i 

to the range space of X then an exact solution exists, but 
is unique only when X is nonsingular in which case it is 
given by the usual inverse X~ X as u = X X v. Otherwise 
there are more linearly independent equations than unknowns 
and an exact solution does not exist. In either case the 

s' ^ 

pseudo inverse solution a v is the unique best solution 

in the sense that, for axiy other possible solution S, it can 

A , 

be shown [48 3 that ||xu-v||'8(|xu-v|| . Finally, 
there are two other properties of the pseudoinverse that 
will prove useful, both of which are easily proven from the 
basic equations (A. 4) and (A. 5). The properties are 
( XX “ X ) T = XX ~ X and XX ~ X X = X. Together they imply that 

— T T -T — T 

( XX ) ( XX ) = XX . 



(A. 6) 
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A. 4 THE STANDARD Z -TRANSFORM 

A thorough treatment of the z-transform is found in 
Jury [49], and its use in the design of digital filters is 
discussed by Goid and Rader 150']., The z-transform is a 
linear map from a vector space of sequences, f = (f m }, into 
a vector space of functions F(z) of the complex variable z* 
It is defined by 


F(z) - f k z k r (A. 7) 

k=0 


-1 -1 

which is a power series in z , thus z is taken within the 
radius of convergence of the power series. Multiplying 
F(z) by z" 1 has the effect of shifting the corresponding 
sequence {f^} by one index as is seen from 


z" 1 F (z) 


f. z 
k 


-(k+1) 


f , z 
n-1 


-n 


k=0 


n=l 


-1 

Consequently, z is known as the unit delay operator . 

The inverse z-transform can be obtained by multiplying both 
sides of (A. 7) by z m L , integrating around a closed contour 
that encloses the origin of the z-plane and applying the 
Cauchy integral theorem to obtain 



00 


E 

k=0 


■fof , ■"**»*. . 
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When the sequence {f^} represents a stable physical 
quantity, for example a sampled voltage waveform, then it 
must not diverge. If F(z) is rational and contains terms 
of the form z/(z- a ), then stability requires that |a| < 1, 
or equivalently, that all poles of F(z) lie within the unit 
circle of the z-plane. Since this result can be extended 
to all stable sequences, the unit circle is the smallest 
contour that always encloses all singularities of F(z); 
consequently, the unit circle will always be taken as the 
path of integration in the z-plane. 

There are two important properties which will be 
used extensively herein; namely, the multiplication 
relation and the Complex Convolution Theorem , known also 
as Parseval’s Theorem. Multiplication in the z-transform 
domain is given by 

00 03 

H(z) = F(z) G(z) = £ £ f n g m z -(m+n) 

m=0 n=0 


oo k 


- Z Z 


f g 
n y kr,n 



so that from (A. 7) it follows that 


k 



n=0 



(A. 8) 


This last expression is the discrete convolution of the 
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sequences {f^} and {g^}. The Complex Convolution Theorem 
can be obtained by considering the sequence ’ (d m > obtained 
as term-by-term products of the two sequences {e m > and 
{f }, i.e., d k ^ e j c f k* The z-transform of {d^} is 


D(z) - 


e k f k z 


-k 


e k z 


-k 


k-1 


skf F(v) v “ * dv 


k-0 


k=0 


= 2njf 


F (v) 


£ e k (z/v) 


-k 


k=0 


v -1 dv 


= ^ij^F(v) E(z/v)v 1 dv, (A. 9) 

which is known as the complex convolution of E(z) and F(z). 
A particularly useful special case is obtained by first 
letting z - 1 to obtain 

oo 

E e k f k = mf FM E(1/v)v " ldv ' 

k=0 

then setting = e k for each k, and finally replacing the 
dummy variable v with z. The result is 


E(z -1 )z“ 1 dz. (A. 10) 

k=0 


This equation reduces to the classical form of Parseval's 
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Identity when the path of integration is the unit circle. 
This corresponds to letting z^ej 10 ' which gives 

00 7T 

E e k 2 = 2¥:. / * ( 3 ®) *(-!“)*• r 

k=0 -TT 

where E(jcu) is the Fourier transform of the continuous 
function e(t) from which {e^} is obtained by sampling at 
the normalized rate of 1 sample per second. Equivalently 
w can be considered as an angle defined by w = ftT, where 
Q is the actual radian frequency and 1/T is the sampling 
frequency . 



APPENDIX 5 


SUPPLEMENTARY RESULTS FOR CHAPTERS 2 AND 3 

A number of statements from Chapters 2 and 3 are 
proven in this appendix. Although all of the proofs are 
important, inclusion in the chapters would have distracted 
from the chapter continuity. Most results are given in 
the "statement-proof ” form, and the section from which the 
statement was taken is found in parenthesis after the proof 
number, e.g., "B.l (2.2)" indicates that the first proof in 
Appendix B is of a statement in Section 2.2. 

B.l (2.2) The matrix A defined by 

a = Jjr $ X ( 2 ) F_(z) F n (z -1 ) Xtz" 1 ! z ' 1 dz 

mn 2Tfj j m n 

is hermitian and positive semidefinitei. 

Proof: By letting z * it follows that 

7T 

a mn " 2ijJ *(J“>X("j“)F m (j»)P n (-ju)j do) 

-TT 

TT 

m hr J I x( i“>| 2 V 3 “ )p n* (j “ >d “ = a nln * ( 

“IT 
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which implies that A is hermitian. To show positive semi- 
definiteness let c be an arbitrary M-vector, then 

c T A c = | c T V(z) V T (z -1 ) c z" 1 dz. 

Since Y(z) = c T V(z), it follows that 


c T A c = *kr f Y(z) Y(z“ 1 ) z " 1 dz 


£ A c = 


IT 


= / I Y <i“> I 2 du i o 


-7T 


Q.E.D. 


B.2 (2.3) The matrix A defined by taking F_^(z) = z m in 
(2.5) is positive definite. 

Proof: Referring to Fig. 2.2, the filter output can be 

written as 

M-l i M-l 

E c m*" m i x(z > 

m=0 


Y(z) = E C m P m (z) XCzl " 

m=0 


= C ( z) X ( z) , 


Where C(z) is defined to be the term in brackets. Using 
this notation the quadratic form c T A c can be written 


„ 1 
S. = 2tF3" , 

f Y(z) Y(z -1 ) z* 1 dz 


IT 

1 

~ 5 ? 

/ 1 C(j<a) | 2 | x (jo) | 2 du> 


-ir 
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Now C(jw) is a linear combination of M complex exponentials f 
and hence must always be nonzero on any finite interval 
tinless c = furthermore, the input signal is converted to 
digital form by an A/D converter operating at a sampling 
rate which is at least twice the highest signal frequency 
in accordance with the uniform sampling theorem. Conse- 
quently, the spectrum X(jo>) must be nonzero on some 
interval of the primary region [-7r,ir], and therefore , 
c T A c > 0 for c 4 £ and A is positive definite. Q.E.D. 

B.3 (2.3) One-step convergence can be obtained for the 
transversal equalizer using the more general algorithm 

c [k+l] = c [k] - 1/2 r[k] VJ (c (k] ) 

where r[k) = diag(Y Q [k3, W > • • • M ) . 

Proofs A derivation analogous to that of (2.14) for this 
new algorithm gives 


k-1 

elk] = J“f .11 - rti}. Aj eior •" 
i=0 

Letting the definitions of A and Q remain unchanged, the 
result corresponding to (2.15) is 

M-i _ „ k-1 

J cor I* = E ^ m L )2 % TT U - Y m Ii] V?, 

m=0 i=0 
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which reduces to zero on the first iteration upon adjusting 

the diagonal components of £[0] to satisfy Y m CO] = 

for m » 0,1, ...M-l. „ 


B.4 (2. 4)' For equation (2.20), the vector b belongs to the 
range space R (A) of A. 

Proof: Since b = F T b and A = F T A F, it is sufficient to 

show R (A) « R(F T ) . Since A is nonsingular and real 
symmetric it has a square root A^ 2 which is also nonsingular 
and real symmetric? consequently, A can be written as 
A * (i 1 /*,)* (A^ 2 F) = S T S. Also, b can be written 

b a F T (A 1//2 a" 1/2 ) b a S T (A -1 ^ 2 b) . Upon observing that 
S T is the adjoint of S it follows from A.l that R(S ) * 

R(S T S) ==' R(A) . However, b is obtained by applying 
S T a f T A^ 2 to the arbitrary vector A ^ 2 b and since 
A 1 / 2 is one-to-one, then R(S T > a R(F T ). Consequently, 

R(F T ) = R (S T ) = R (A) . q.e.D. 


B.5 (2.5) The map P = XX*" 1 defined by (2.33) is a projection 
matrix, i.e., it is both self-adjoint and idempotent. 

Proof: First, P is hermitian since it is real and 


T 

P = 


N-l 

E 


N-l 


= E C = ^ • 


m=0 


m=0 
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Furthermore, since 


p 2 = 


N-l 

\ 

/ N-l 

E 

4^ 1 

E *n ^ 

m=0 

/ 

\ n=0 

N-l 

N-l 


E 

T in 

-%i 

T 

6 i = P 

mn x n — r 

m=0 

n=0 



then P is also indempotent, hence P represents a projection 
operator. Q.E.D. 

B . € (2.5) The matrix £ T P ¥ - diag.tl^ | 0^^) * 

Proof: The result follows directly from the fact that {^} 

is an orthonormal set: 


p f - 


to 

h 


T 




T 


N-l 

E In ^ [V 4.i* • * • 4*i-i [] 

m=0 


& 


= diagCl^ | ^ ) . 


— M-N 


Q.E.D. 
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B.7 (3.1) If a priori statistical information is available 
for evaluation of the gradient (3.4), then the tap-gain 
adjustment algorithm (3.2) reduces to other well-established 
algorithms . 

In the first case it is assumed that the noise 
statistics are unknown, the transmitted signal and noise 
are uncorrelated, and that statistical information is known 
about the transmitted sequence {d[k] } . Under the decision- 

A 

directed hypothesis that d[k] = dlk] for all k, the 
gradient (3.4) becomes 

^ VJ(cIk]) = B{v [k] v T [k]} elk] - f 

where ct^3 i- 6 assumed independent of v[k] and 

r dv - r <jv (n) A E{v[k] dlk]} - P T E{xlk] dlk]}. 

The result, known as the Griffiths Algorithm, was considered 
by Griffiths [51 ] in connection with antenna-array signal 
processing, and some convergence results were developed. 

For the second case it is assumed that complete 
statistical information is available a priori, so that (3.4) 
becomes 

| VJ(c[k]) = R vv £[k] - r dv , 

with r, defined as before and the M x M correlation matrix 
— dv 

R vv given by R yv = R VV (Q) = E{v[k] v T [k]}. The result is 
the well-known steepest descent algorithm 


c[k+l] = c Ik] - Y (Syy c Ik] - rj v ) . 
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B.8 (3.2) Under the stated hypothesis the matrix R xx is 
positive definite. 

Proof: Since xlk] « s [k] + ntk] for all k, then the 

signal and noise vectors s_[k] and n [k] , respectively, 

corresponding to xjk] are defined component -wise by 

xlk] = J3 Ik] + njkj . It follows from the independence of 

signal and noise that R = E {sjk] Ik ] } + E {n[k] n T [k] } 

^ R + R Now R is a correlation matrix, hence positive 
— *s s — nn — s s 

semidef inite ; also R nn can be diagonalized by a unitary 

T 

matrix Q to give Q R xx Q = A . However, upon defining 
A T 

m « q n a new representation of the noise results with 

T T T 

crosscorrelation matrix E{Q n[k]n IK] Q > - Q R nn Q = A^, 

i.e., the crosscorrelation matrix of m[k] is diagonal. , 

The elements are all variances, thus are all positive, 

and it follows that R and R are positive definite. 

— nn —xx L 


Q.E.D. 



APPENDIX C 


COMPUTER SIMULATION 

This appendix contains a detailed description of 
the digital computer simulation written in FORTRAN to 
study the equalizer design which has been developed. 

The simulation consists of a main routine, called MAIN, 
and a number of subroutines, which together simulate a 
transmitter, digital communication channel, and the 
equalizer. In addition to using subroutine GRAPH for 
producing graphical output of the equalizer performance, 
MAIN uses two other subroutines, namely FLTSET and HADGEN, 
to compute the initial tap-gain settings in accordance 
with (2.37). Each of these in turn calls two or more 
subroutines as shown in Fig. C.l. MAIN is described in 
C.l while the remaining subroutines are covered in C.2. 

Two numerical methods for finding eigenvalues and eigen- 
vectors of real symmetric matrices are described in C.3. 

One is used by the IBM library subroutine EIGEN (52 J and 
the other by the UNIVAC library subroutines TRIDMX, EIGVAL, 
and EIGVEC (53], It is interesting to note here that, 
since a modern digital computer is actually a very 
sophiscated digital filter in itself, then the simulation 
of the equalizer is quite realistic, although not very 
efficient from a hardware standpoint. 

1CL0L 
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C.l PROGRAM MAIN 

The program MAIN can be subdivided into two basic 
parts. The first simulates the transmission of an 
initializing pulse through the channel and performs the 
equalizer initialization based on the theory of Chapter 2. 

The second part generates the transmitter message sequence 
and uses the initialized equalizer to track the channel. 

A flowchart for MAIN is shown in Fig. C.2, and the program 
listing in Fig. C.3. For reasons of continuity a few 
parameters , usually pertaining to small blocks of code, 
are set locally as they occur rather than at the beginning 
of the program. The listing of key parameters found in 
Table C.l will be helpful in identifying and using these 
variables. 

A concise description of the program is found in 
block 0. The program is designed to repeat until no 
further data cards are found, at which time it terminates 
normally. After reading the simulation parameters and 
initializing variables in blocks 0-30, the channel output 
sequence (X) is formed in blocks 40 and 50 by passing an 
initializing pulse through the noisy channel (CHNL, G, 
CNOISE) . Block 60 calls the subroutine FLTSET which returns 
in common the eigenvalue and eigenvector information (FSET) 
for computing £2.37). In block 70 an M x M Hadamard matrix 
(H) is formed by subroutine HADGEN, and the matrix (F) , whose 
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Fig. C.2. Simulation Program MAIN 














Fig. C.2 


(Continued) 
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DIMENSION U< 32) *Xt 52 ) *G(^2) *CHNL<32> * IKEG< 10) r C t 32? ) r 
1 iw ( 32 ) * R ( 32 ) * F ( 32 ► 32 ) * UELC ( 32 ) > TRANS ( 32 > * PLOT Y ( 64 ) 
INTEGER DELAY *UMAX*TRTYPE *TOTLER 
COMMON H<32*32> *FSET(32*32> 

LOU I VALENCE ( DLLC ( 1 ) * CHNL ( 1 ) *H(1* 1) ) » < TRANS (1 ) »H(1*2> ) 
1 * { PLOTY (1) *H ( 1 * 3) ) 

C 

COOO THIS PROGRAM SIMULATES AN ADAPTIVE DIGITAL EQUALIZER 
C AND ASSOCIATED BASEBAND CHANNEL. THE INITI AL I Z AT I ON 

C ‘ IS PERFORMED IN BLOCKS 10-90* AND ON-LINE OPERATION 
C IS SIMULATED BY BLOCKS 100-150. 

C 

c 

c THE DATA CARD FORMATS ARE AS FOLLOWS 

CARDl M < 1 5 ) * N(I5)* DELAY ( 15) * MSGLTH(I5>* TRTYPE ( A6) 

CARD2 CMEANIF10.5) * CSTUEV (F10.5) * CVFCTR (FlO . 5 ) * NINITCI5) 
CARD3 CHANNEL IMPULSE RESPONSE LENGTH NG(I5> 

CARD4 CHANNEL IMPULSE RESPONSE G(8F10,5) 

C 

C REPEAT ENTIRE SIMULATION UNTIL DATA RUNS OUT 

DO 500 11=1*15 

READ (5* 1 *Ei jO=5U1 ) M*N*DELAY*MSGLTH* TRTYPE * 

1 CMEAN* CSTUEV r CVFCTR * NIN IT 
1 FORMAT ( 4 i 5 * 4X * A6/ 3F 10 * 5 * I 5 ) 

C 

C INITIALIZE VARIABLES AND ARRAYS 

DO 10 1=1*32 
CHNL ( I ) =0 
C(I)=U 
D(I)=U 
G(1)=U 
XU)=0 

DO 10 J=1 * 32 
F ( I * J ) =0 
H( I* J)=0 
10 CONTINUE 

C SET THE NOISE GENERATOR PARAMETERS 

IA=127329 
IB=566387 
IC=2**20 
REALC=IC 

CVAR=CSTUEV*CSTDEV 

C PROGRAMMING NOTE: CSTDEV MUST NOT BE ZERO 

SNR= 20*ALOG10(1/CSTDEV) 

C PRINT THE INITIAL EQUALIZER PARAMETERS 

WRITE (6*12) M* Mr MSGLTH* TRTYPE * CMEAN * C VAR * SNR * DEL AY * 

Fig. C,3. Listing of Program MAIN 
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c 

cozo 

c 

21 

22 

C 

C030 

C 

C 

C 


c 

CU40 

C 

C 


42 

C 

41 

C 

CUSO 

C 

5U 

52 


53 


Fig. C 


1 LVFCTKfNlNIT 

F-UKMAT (1H1////10X# ’SYSTEM PARAMETERS’/ 

1 1UX » ’NUMBER OF FILTER SECTIONS (M> = f »l5/ 

2 tOX » * LENGTH OF FILTER SECTIONS (N) =*»I5/ 

3 10X f * LENGTH OF TRANSMITTED MESSAGE (MSGLTH) t 17 / 

4 1 UX » ’ TYPE OF TRANSMITTER SIMULATED (NTRTYP) =’*Ao/ 

5 10X # ’ NOISE MEAN = * t F 10 * 5/ 1 OX r ’ NO I SE VARIANCE =**F10 

6 1 OX r ’ S I GNAL TO NOISE RATIO IN DB =*»F10.5/ 

7 10X » * TOTAL DELAY = ’ » lb / 10X »• CONVERGENCE FACTOR = 
B F 10 .5/iUXr ’NUMBER OF INITIALIZING PULSES =»»Ib) 

READ IN THE CHANNEL DESCRIPTION 

READ (5*21) NO f (G( I ) t 1 = 1 *NG) 

FORMAT ( 1 5/ ( 8F 10 • 5) ) 

WRITE ( ta r 22 ) NG * ( I » G ( 1 ) * I = 1 * NG ) 

FORMAT (////10X* ’CHANNEL RESPONSE G OF LENGTH’ » IS/ 

1 (10X* ’G( ’ * 12* ’ ) = » *F8.4) ) 


IF INITIALIZING SEQUENCES OTHER THAN THE UNIT-PULSE 
ARE USED THEY SHOULD bE DEFINED IN THIS BLOCK* 
OTHERWISE A UNIT-PULSE SEQUENCE IS GENERATED 

rju=i 

U ( 1 ) -1 

ADJUST THE INPUT SEQUENCE LENGTH 

CHECK FOR EXCESSIVE CHANNEL-RESPONSE LENGTH 
L-NO+NG-1 

IF (L .LE. M-N+l) GO TO 41 
WRITE (bf 42) 

FORMAT ( 1H1 » 10X* ’ ***THE OUTPUT IS TOO LONG***’) 

GO TO S00 

EXTEND THE OUTPUT SEQUENCE LENGTH 
L-M-N+l 

FORM THE CHANNEL OUTPUT SEQUENCE X ( Z ) =D ( Z ) *G ( Z ) 

CONTINUE 

DO S2 1=1 *NG 

CHNL l 1 ) =U 

DO 51 K = 1*L 

NGM=NG-1 

DO 53 l=i#NGM 

CHNL(NG-I+1)=CHNL(NG-I ) 

CHNLU )=D(K ) 

.3. (Continues! Lifting of program ftAiN 
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b 1 * 

C 


bb 

C 

bl 

C 

bo 

C 

COoO 


070 


7 U 

c 

CUtiO 

c 


ai 

aO 

c 

CU90 

C 


Pig. 


X(K)=0 

UO b4 I — i f NO 
X(K)=X(K)+G( I ) t-CMNLU) 

GENERATE A GAUSSIAN RANDOM VARIABLE FOR CHANNEL NOISE 
KSUM=0 

DO 5b 1=1,12 

Ib=MOD( lB*(2**10+3> » IC) 

REALB=10 

RSUM=RSUM+REALb/REALC 

CN0ISE=(RSUM-6)*CSTDEV/NINIT + CMEAN 
ADD THE CHANNEL NOISE 
X(K)=X(K)+CNOlbE 
CONTINUE 

PRINT THE CHANNEL OUTPUT FOR THE UNIT-PULSE 
WRITE (6rS6) Lr ( I » X ( I ) , 1 = 1 »D 

FORMAT (////10X» ’CHANNEL OUTPUT SEQUENCE X OF LENGTH ’ t 
1 I5/(10X’X( * , 12, ’ )=’ *F8.4> ) 

COMPUTE THE FILTER SETUP MATRIX FSET 

SELECT TYPE OF EIGENVALUE COMPUTATION 
EIGTYP=’ JACOBI* 

CALL FLTSET ( M * N , X t E IGTYP ) 

GENERATE AN MXM HADAMARD MATRIX AND NORMALIZE 

CALL HADGEN(M) 

RLALM=M 

SQM=bQRT (REALM) 

DO 70 1=1 »M 
DO 70 J=1 » M 
Hi I» J)=H( I » JJ/SQM 
CONTINUE 

SET THE FILTER SECTIONS 

DO 60 J= 1 , M 
DO 80 1 = 1, N 
F ( I , J ) =0 
DO 81 K = 1,N 

F(I»J)=F( I, J)+FSET(I»K)*H(Kf J> 

CONTINUE 

INITIALIZE THE FILTER AND SET THE TAP GAINS 

DO 91 1 = 1 * M 
DELC ( I ) =0 


C P 3. (Continued) Listing of Program MAIN 
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91 kll)=0 
AMSE=0 
UO 9 id K - 1 * M 

C LOAD AND SHIFT THE FILTER REGISTERS 

NM=N-i 

DO 93 1=1 rWM 

96 K(N-1+1 )=« (N-l ) 

R ( 1 ) =X ( K ) 

0 FORM THE FILTER OUTPUT Y 

Y = 0 

DO 94 I = 1 * M 

W(I)=0 

DO 96 J=1 9 N 

9b W(I)=W(I)+K( J)*F(J*I) 

94 Y=Y+vm>*C(I> 

C COMPUTE THE BIT ERROR AND THE TAP ADJUSTMENTS 

E = Y 

IF (K *GT .DELAY) E=Y-D ( K-DELAY ) 

AMSE=AMSE+ < E*E- AMSE ) /K 
DO 96 I=1*M 

9b DELC ( I ) =DELC (I ) +E* W ( I ) 

9 id CONTINUE 

L SET THE EQUALIZER TAPS 

DO 97 I = 1 * M 

97 o(l')=C< I )-DELC ( I ) 

WHITE (b*98) ( I r C ( I ) * I = 1 * M ) 

9b FORMAT ( ////iOX * ’INITIAL TAP GAIN SETTINGS * / 

1 (10X* *C ( » f 12* * > = ' *F8.4) ) 

WRITE 16*99) MSGLTH 

99 FORMAT ( 1,11 * 25X * * THE TRANSMITTED SEQUENCE IS»*I4* 

1 » BI IS LONG * ////6X * 11HTRANS/DEC ’ N* 9X * ’NOISE* * 4X * 

2 *CHNL OUT* * 2X * * EQLR OUT ♦ * 4X * * ERROR * * 3X * * A VE MSE » * 

3 2X*' TOTAL MSE * * 4X * * C ( I ) >*//) 

luO THE REMAINING PART OF THIS PROGRAM IS A LARGE LOOP 
THAT ITERATES ONCE FOR EACH MESSAGE BIT 

110 SHIFT TO RUN MODE AND TRANSMIT THE MESSAGE . 

DO 111 I = 1*L 
TRANS ( 1 >=0 

111 CHNL ( I ) =0 

DO 113 1=1*64 
113 PLOTY(l)=0 
AMSE=0 
TMSE=0 

totlek=u 


Fig* Co 3, (.Continued) Listing of Program MATN 
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NC OUT-4 
LL-8 
LLM=LL-1 
UMAX=25 
ELIM=2 * 0 
NGM=NG-1 
NLINt-1 
C 

DO 112 K=1 * MSGLTH 
C 

C 120 GENERATE A TRANSMITTER PULSE 

C SELECT THE INPUT SEQUENCE FROM EITHER A PN 

C SEQUENCE (PNSEQ) OR RANDOM NUMBER GENERATOR (RANDNO) 

IF (TRT YPE • EQ » ■ PNSEQ * ) GO TO 122 
IF (TRTYPE.EQ. ’RANDNO 1 ) GO TO 121 
C ASSUME THAT A UNIT-PULSE IS DESIRED 

DD-0 

IF (K.EQ-1) DD=1 
GO TO 126 

121 CONTINUE 

C GENERATE A RANDOM NUMBER BETWEEN 0 AND 1 USING 

C THE POWER RESIDUE METHOD 

1A=MG0 (IAM2**10 + 3) # IC) 

REALA-I A 

RANDIN=REALA/REALC 

C CONVERT TO BINARY USING A THRESHOLD OF 0.5 

DU=0 

IF ( RAND IN .LT. G.5> DD=1 
GO TO 126 

C GENERATE A 2**LL-1 BIT PN SEQUENCE FOR THE MESSAGE 

122 IF (K.NE.l) GO TO 123 
DO 125 I=i*LL 

125 I REG ( I ) =1 

123 DD=IREG(3) 

ITEMP=MGD< IREG(2) + 1 REG ( 3 ) + 1 REG ( 4 ) +IR£G<8) »2) 

DO 124 1=1 r ELM 

124 ireg(ll-h-d = ireg(ll-i) 

IREG ( 1 > = I TEMP 

12b CONTINUE 

STORE THE TRANSMITTED SEQUENCE FOR PRINTOUT LATER 
DO 127 1 = 1 f DM AX 

j.27 TRANS (DMAX-I + 1)=TRANS(DMAX-I ) 

TRANS ( 1 ) =OD 

130 PASS THE TRANSMITTED PULSE THROUGH THE TIME-VARYING 


Fig, C>3, (Continued) Listing of Program HA,IN 
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C CHANNEL 

c 

100 CONTINUE 

C THE TIME-VARYING CHANNEL SHOULD BE CHANGED HERE 

DO 101 1 = 1 p N6M 

101 CHNL(NO-I+l )=CHNL(NG-I ) 

CHNL ( 1 ) =DU 

XX = 0 

Do 102 1 = 1 p NG 

102 XX=XX+G< I ) *CHNL< I ) 

C 

C GENERATE A GAUSSIAN RANDOM VARIABLE WITH PRESET 

C MEAN AND STANDARD DEVIATION 

RSUM=0 

DO 100 1=1 p 12 

lb=MOD( IB*(2**10+0) »IC) 

RLALb=IB 

100 KbUM=KSUM+KEALB/REALC 

CHO I SE= ( RSUM-6 ) *CS TDEV+CMEAN 
C ADD THE CHANNEL NOISE 

XX=XX+CNOISE 
C 

C140 FILTER THE CHANNEL OUTPUT PULSE XX 

C LOAD AND SHIFT THE FILTER REGISTERS 

DO 141 1=1 pNM 
141 R(N-I+1>=R<N-I> 

R ( 1 ) = X X 

C FORM THE FILTER OUTPUT Y 

Y=0 

DO 142 1=1 pM 
W Cl )=0 

DO 140 J=1pN 

mo W ( I )=W ( I ) +R ( j) *F ( Jp I ) 
m2 Y=Y + W( I )*C< 1) 

C SET UP THE EQUALIZER OUTPUT FOR PLOTTING 

I=2Q*( Y+0.5>+l.B 

IF(I.LT.l) 1=1 
IFU.6T.S1> I = S1 
PLOTY(I)=PLOTYU)+l 
C 

Cl SO COMPUTE THE ON-LINE TAP-GAIN ADJUSTMENTS 
C 

C ESTIMATE THE TRANSMITTED VALUE 

ESTU=0 

IF(Y.GE.O.b) ESTD=1 
E=Y-ESTD 


fig. C . 3 . (Continued) Listing of Program JW1N 
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AM5E=AM5L+ (E+E-AMbE) /K 
TM5t=TMSE U. + E 

C CHECK FOR U I VLRGENCE 

IF ( ADS ( E. ) .LT. LLIM) GO TO 153 
WRITE ( 6 * 151) K 

151 FORMAT (////10X* »***THE ALGORITHM DIVERGES AT STEP* * 
1 1 5 » f ♦♦♦ * ) 

GO TO 500 

C ADJUST THE FILTER TAPS USING THE GRADIENT ESTIMATES 

153 CO 152 I = i*M 

152 C ( I ) -C ( I ) — E*W ( 1 ) *C VFCTR 

C PRINT EVERY NLINE-TH LINE FOR BREVITY 

IF (K.NE.K/NLINE+NLINE) GO TO 155 
WRITE (6*154) K * TRANS ( 1+DELAY ) ► ESTD * CNOISE * XX * Y *E * 

1 AMSE * TMSE* (C(I) *I=l*NCOUT) 

154 FORMAT ( 1H *I6*F4.i** / * * F4 . 1 * 6X * 4F10 . 5* 2610 . 4 * * *•* 
1 4F10.5/(64X*4F10.5> ) 

155 CONTINUE 

IF ( ABS (ESTD-TRANS ( 1+DELAY ) ) .LT. 0.1) GO TO 157 

TQTLER=TOTLER+l 

WRITE (bflbu) 

lbb FORMAT (lH+rl X*»*+*** f ) 

157 CONTINUE 
C 

112 CONTINUE 
C 

0500 TERMINATING PROCEDURE 
C 

WRITE (6*502) TOILER * MSGLTH 

502 FORMAT 1//26X * * THERE WERE** 14* * ERRORS IN THE* *14* 

1 • BIT MESSAGE*) 

CALL GRAPH ( -0 . 5 * 0 . 05* PLOT Y *51) 

500 CONTINUE 

501 CONTINUE 
STOP 
END 


Fig. C.3. {.Continued} Listing of Program MAIN 
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TABLE C . 1 

LIST OF SIMULATION VARIABLES 


Name 

T2E£ + 

Description 

AMSE 

R 

Average MSE per iteration 

CHNLC32) 

R 

Shift-register for simulating the channel 

CMEAN . 

R 

Mean value of channel noise 

CNOISE 

R 

Channel noise value 

CSTDEV 

R 

Standard deviation of channel noise 

CVAR 

R 

Variance of channel noise 

CVFCTR 

R 

Convergence factor for the on-line 
algorithm 

D(32) 

R 

Working array for the initializing 
sequence 

DELAY 

I 

Delay of reference pulse in bits 

DELC (32 ) 

R 

Tap-gain increments for initialization 

DMAX 

I 

Maximum delay 

E 

R 

Equalizer output error 

ELIM 

R 

Equalizer output error threshold 

ESTD 

R 

Decision on equalizer output 

F ( 3 2 , 3 2 ) 

R 

Matrix whose columns are the filter- 
section tap-gains 

FSET C32 , 32 ) 

R 

Matrix returned by FLTSET used for 
computing initial filter-section 
tap-gains 

G (32) 

R 

Noise-free channel unit-pulse response 

H (32 , 32 ) 

R 

Hadamard matrix returned by HADGEN 


•j* 

R indicates a real variable and I indicates an 

integer. 
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Name 


TABLE C . 1 (Continued) 


area 


Description 


L 


I Channel output sequence length 


M 


I Number of filter sections 


MSGLTH 


I Number of bits in transmitted message 


N I 

NCOUT I 

NG X 

NINIT I 

NLINE I 

PLOTY (60) R 

R (32 ) R 

SNR R 

TMSE R 

TOTLER I 

TRANS (32 ) R 


Length of shift register 

Number of tap-gain values to be 
printed at each iteration of the 
on-line algorithm 

Length of noise-free channel unit- 
pulse response 

Number of isolated initializing 
pulses used to set the equalizer 

Number of output lines skipped before 
printing a line of output data 

Contains the equalizer output-value 
distribution for plotting by GRAPH 

Shift-register for equalizer filter- 
sections 

Channel SNR on a per-bit basis 

Total (cumulative) MSE 

Total number of errors in equalizer 
output 

Array to store transmitted sequence 
for later printout 


TRTYPE Transmitter type, either ’ PNSEQ' or 

' RANDNO ' for PN sequence or' random 
number generator, respectively ? otherwise 
a unit-pulse sequence is transmitted 


W (32) 


R Filter-section output vector 


X (32) 


R Channel output sequence 


Y 


R Equalizer output 


V 
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columns are the filter-section tap-gains , is then computed 
as the product of FSET and H. The equalizer is now 
initialized in block 80 by first computing the equalizer 
tap-gain corrections (DELC) , and then adding the result 
to then existing tap-rgains (C) in accordance with (2.2), 
and the initial tap-gains are printed in block 90. 

The remaining part of MAIN, starting with a 
description in block 100, is essentially a large loop 
that iterates once for each transmitted pulse. All coding 
inside the loop has been done inline, i.e., there are no 
subroutines, in order to reduce the program execution time. 
After a number of variables are initialized in block 110, 
the loop begins. A transmitter pulse (DD) is generated 
in block 120 using either a pseudo-noise sequence (TRTYPE « 
‘FNSEQ 1 ) or a random binary number (TRTYPE = 'RANDNO') gen- 
erated by quantizing random numbers that are uniformly dis- 
tributed between 0 and 1. If no transmitter type is specified 

then a unit-pulse is transmitted. - In block 130 the pulse 
is next passed through the channel which is corrupted by 
gaussian noise. Since time -variation in communication 
channels is rather arbitrary, this effect is simulated by 
inserting a user-supplied description immediately after 
statement 130. The transmitted pulse (XX) is then filtered 
by the equalizer in block 140 to produce the equalizer 
output (Y) . The average MSE per bit (AMSE) and total 
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(cumulative) MSE (TMSE) are next calculated, and the bit 
error (E) is compared against the maximum acceptable error 
value (ELIM) . If the limit is exceeded an error message 
is printed and the equalizer can either be reinitialized 
by returning to block 50, or the run can be terminated 
by a transfer to 500. Otherwise, at block 150 the 
equalizer tap-gains are adjusted in accordance with (3.6) 
and, when indicated by a local parameter (NLINE) , a line 
of data concerning the iteration is printed before 
repeating the loop for the next transmitter bit. After 
the indicated number (MSGLTH) of message bits has been 
transmitted, the total number of bit errors is printed 
and GRAPH is called to plot the frequency of occurrence 
of each value of the equalizer output. 

C . 2 SUBROUTINES 

The flow chart and listing for subroutine FLTSET 
are found in Figs. C.4 and C.5, respectively. It generates 
the convolution matrix (XMX) of (2.2B) from the channel 
unit-pulse response vector (X) that is passed from MAIN. 

The N x N real symmetric matrix (SYMX) is then formed and 

/ 

the corresponding eigenvalues and eigenvectors determined 
by one of the two methods discussed in C.3. If the Jacobi 
method is used (EIGTYP = 'JACOBI 1 )# then subroutine MODE 
is used to convert SYMX into the linear array (EVAL) form 
required by the IBM library subroutine EIGEN. The other 
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Fig. C.4. Subroutine FLTSET 
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SUbROUTlNc FLTSET ( M > N » X > L I G T Y P ) 

DIMENSION X (32) » XMX (32*32 ) > EGVL < 32 ) » D < 32 ) > SD ( 32 ) > 

X TEMPI ( 32 ) > TLMP2 ( 32 ) * SYMX ( 32 > 32 ) * EVAL ( 528 ) » EVC7 ( 1 U24 ) 
COMMON EGVC (32 >32) *FSET (32 >32) 

EQUIVALENCE (SYMX(1>1) >FStT(l*l) >EVCT(1) ) 

EQUIVALENCE (EVAL(l) >EGVL (1) ) > (EVAL(33) >D(1) ) » 

1 ( EVAL ( 6S ) > SD ( 1 ) ) » ( EVAL ( 97 ) *TEMPU( I) ) > (EVAL( 129) > 

2 TEMP 2 ( 1 ) ) 

FORM THE CHANNEL OUTPUT MATRIX XMX 

DO 1 1=1>M 

DO 1 J- 1 # N 
XMX(I>J)=0 

IP ( X • GE • J ) XMX ( I » J ) -X ( I-J+ 1 ) 

CONTINUE 

FORM THE REAL SYMMETRIC MATRIX SYMX= ( XMX ) T*XMX 
DO 2 X=1»N 
DO 2 J=1>N ' 

S Y M X ( X > J ) - 0 
DO 3 K-l * M 

SYMX( I> J)=SYMX(I> J)+XMX(K* I )*XMX(K> J) 

CONTINUE 

IF (ElGTYP.EU. * JACOBI • ) GO TO 10 
IF (EIGTYP.EQ. f C0D1AG* ) GO TO 5 
OTHERWISE PRINT DIAGNOSTIC AND GO TO 10 
WRITE (6>b) 

FORMAT (10X> »****EIGNEVALUE METHOD NOT SPECIFIED***** > 
GO TO 10 

CONTINUE 

COMPUTE THE EIGENVALUES AND EIGENVECTORS OF FSET 
CALL TR IDMX ( N > 32 » SYMX » D > SD ) 

C ALL EIGVAL (N> EGVL*U*S0 > TEMPI * TEMP2 ) 

CALL EIGVEC(N>32*SYMX>D#SD* EGVL t EGVC > TEMPI * TEMP2 ) 

CALCULATE THE FILTER SET-UP MATRIX FSET 
DO 4 J“ 1 # N 

EGVL ( J) =1 /SORT (EGVL( J) ) 

DO 4 i=l»N 

FSETi Ir J)=EGVL( J)*EGVC (I >U) 

4 CONTINUE 
RETURN 
C 

10 CONTINUE 

F±g f C . 5 . Listing of Subroutine FLTSET 
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C 

CALL MODE (li» SYMX » LVAL ) 

CALL LlGCN ( LVAL t LVCT t N • 0 ) 

CALCULATL THt FILTER StT-UP MATRIX 
DO II JJ=1 f M 

If- (J.LE.N) GO TO 12 
DO 13 1=1 fU 

10 FSET ( I » J ) =0 
GO TO 11 
K=(J-1)*N 
KK=J* ( J + l ) /2 

LVAL ( KK ) = i/SON T < EVAL ( KK ) ) 

DO 14 I = 1 * U 

14 FbET( I r J)=EVAL<KK)*EVCT (K+X) 

11 CONTINUE 
RETURN 
END 


. tContinuedX Listing ot Subroutine FLT SET 


Fig, C . 5 
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method (EIGTYP » 'CODIAG* ) is that of Householdei , which 
uses the UNIVAC library subroutines TRIDMX to convert SYMX 
to tridiagonal form, then EIGVAL and. EIGVCT to compute the 
corresponding eigenvalues (EVAL) and eigenvectors (EVCT) . . 
These subroutines require several working matrices (TEMPI >" 
TEMP2 , SD, D) • After the eigenvalues and eigenvectors are 
determined, FLTSET computes the filter set-up matrix (FSET) 

and returns. In the notation of 2.5, the columns of FSET 

—1 ' 
are given by £-k' k = 

Subroutine MODE places the upper triangular part of 
a given real symmetric matrix (SYM of dimension DIM) in a 
linear array (EVAL) columnwise starting with column 1 as 
required by subroutine EIGEN, then returns. The correspond- 
ing listing is given in Fig. C.6. 

For any order (M) which is a multiple of four and 
is not greater than 32, subroutine HADGEN generates a 
Hadamand matrix (H) by one of three procedures, depending 
on the order desired. The flow chart is shown in Fig. 0.7 
and the corresponding listing in Fig. C.8. If M is a 
power of 2, i.e., 4, 8, 16, and 32, then H is formed in 
the manner discussed in A. 2 using subroutine DOUBLE to 
repeatedly make the substitutions 


■— - 


r 1 

1 1 


H* 

1 

H 


and -1 


1-1 


l-i lj 

" —B 


L- — J 



r r 


\ 


120 


g 


C 


i 


^r-a%W l e r0 «PV, 


SUbKuUT INL MuUt. (UiMrSYMtLIN) 

ThlS PROCEDURE PLaCES THE UPPER Tk I ANGULAR PART OF 
SQUARE REAL SYMMETRIC MATRIX SYM IN LINEAR ARRAY LIN 
IN MODE 1 FORM AS REQUIRED BY SUBROUTINE EIGEN 


DIMENSION SYM( 52f3Z ) r LIN ( 528 ) 
INTEGER DIM 
REAL LIN 
L=DIM*(DIM+1 ) /£ 

DO 1 I - 1 r L 
L1N(I)=0 
DO 2 J=1»UIM 
DO 2 I=1>J 
K=J*(J+l)/2 

IF (l.LT.J) K- I + ( J~1 ) * j/2 
L i N ( K > - S Y M ( I » J ) 

CONTINUE 

RLTUKN 

END 


C , 6 


Listing of Subroutine MODE 



UADtJEN 


NO 


IS M 



Fig. C.7. Subroutine HADGEN 
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C 
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4 


3 


4 

1 



SUjKOUT iNL HAUGEN (M> 

THIS PROCEDURE GENERATES A BINARY HADAMAKU MATRIX OF 
OkDER M = 4*K vvllh K <-rtr AND C AN EASILY BE EXTENDED 
TO K>8 


Common hi 32 *32) 

INTEGER RtS(2U) 

IM (M/4 ) *4.NE.M.0k.M.GT.32) GO TO 25 

If M-2o FORM A S 1 ZL 14 MATRIX AND DOUT5LE 
IF (M*NE *2o) GO TO 10 
LJ = 13 

CALL CADRES ( LJ r Lk r RES ) 

Lo=LJ4l 

DO 1 I-1*LJ 

DO l J- I r L J 

If ( J.NE. I ) GO TO 2 

111 I r J ) — U 

GO TO l 

lK(I.NL.l) GO 10 5 
Hi I r J) = l 
Hi Jrl)=l 
GO TO l 
DO 4 K“ 1 r LK 

I F ( J” I • £0 • RES ( K ) ) GO TO 5 

CONTINUE 

tld r J)=-l 

H(J#I)=-1 

CO NT aimUE 

CALL DOUBLE (LJ) 

RLTUKN 


IF M— 12 OR M-20 FORM THE MATRIX DIRECTLY 
IF M-24 FORM A SIZE 12 MATRIX AND DOUBLE 
1U if (M*NL* 12. AND.M* I4E. 20* AND# M*NE. 24) GO TO 20 

LJ-11 

IF ( M • EO • 20 ) LJ-19 
CALL CADRES < LJ » LK r RES ) 

LJ=LJ+l 

DO 11 1 — 1 r L J 

DO 11 J- 1 * L J 

IF ( I • NE. • 1 ) GO 10 13 

H( 1 1 J) = i 

Hi Jr 1 ) = 1 

GO TO 11 

13 IF(J.NE.l) GO TO 15 


Fig. C.8, Listing of Subroutine HADGEN 
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00 10 11 

15 00 lb K-l t LK 

IF ( J-l ,EQ. KES(K> ) GO TO 18 
lb CONTINUE 
H ( 1 # J )=-1 
H(J*I)=1 
GO TO 11 

18 H< I t J)=l 
H(Jrl)=“l 

ii continue: 

1F(M*NE.24) GO TO 19 
CALL UOUbLL (LJ) 

19 CONTINUE 
RETURN 

SINCE M=2**K» FORM A SIZE 2 MATRIX AND DOUBLE 
RECURSIVELY K-l TIMES 
2 U H(l»i)=l 

H(i»2)=l 
H ( 2 r 1 ) = 1 
H(2*2)=“l 
Ll = l 

IMM.EQ. 8) L X -2 
IF (M. EG. lb) LI = 3 
IF ( M.EU. 32 ) LI-4 
L-2 

22 UU 24 I - 1 » L I 

CALL DOUBLE (L) 

24 L-2 *L 

RETURN 

2b WR I TE ( 6 * 26 ) 

26 FORMAT (1H1*****THE VALUE OF M IS INCORRECT****) 
RETURN 
END 


Fig. C.8. (Continued) Listing of Subroutine HADGEN 
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If M is not a power of 2 then the generating pro- 
cedure is slightly more complicated, and requires the 

notion of the Legendre symbol . If the congruence relation 
2 

r s q (mod p) has a solution for integers p, q, and r, then 
q is called a quadratic residue of p, otherwise it is 
called a quadratic nonresidue. The Legendre symbol (q/p) 
is defined for a prime p to be 1 when q is a quadratic 
residue of p, and -1 when q is a quadratic nonresidue. 

The remaining values of M now fall into two categories. 

If M * 2 k Cp + 1), k = 0, 1,... , for some prime p = 3 
(mod 4) , i.e., M « 12, 20, 24, then the (p + 1) x (p + 1) 
array H 1 is first formed, with each element computed 
according to 

if i ■ 1 or j ■ 1 

if i r j and i > 1 

otherwise 


h . , = 


f i - 
Cj-i/p) 


where the quadratic residues are determined by QUADRES. 

If M = p + 1, then H is taken to be H* , otherwise 

M = 2 (p + 1) and DOUBLE is called to generate H before 

k 

HADGEN returns. Finally, if M = 2 (p + 1) for prime 
p = 1 (mod 4), i.e., M = 28, then H’ is computed according 


to 


/ 





if i = 1 or j = 1 but i / j 

if i r j 
otherwise , 


1 

0 

(j-i/p) 
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and again QADRES is called to obtain the quadratic 
residues. Since M = 2 (p + 1) , then DOUBLE is called to 
obtain H before returning. A detailed treatment of the 
generation of Hadamard matrices for all order up to 100, 
and many orders greater than 100, can be found in the 
literature (44,45)* Consequently, HADGEN could be expanded 
to handle considerably larger orders than the present 
value of 32. 

Subroutine DOUBLE doubles the size of a given matrix 


H' = Ih’*,] with elements 1, -1, and 0, to obtain a new 

ij 

matrix by making the following substitutions. 



’l 1 


r — — 1 

1 -1 


i 1 
*-* 

i 

h ij = ^ 


, = 0* 


, and h|j - -!♦ 



1-1 


-1 - 1 _ 


-1 1 


The original matrix is overwritten by the new matrix H 
which is binary, i.e., all elements are 1 or -1. The 
listing is found in C.9, and the algorithm is to simply 
test each element h^ and make the corresponding matrix 
substitution just indicated. 

Subroutine QADRES computes the quadratic residues 
of a given prime defined in the discussion of HADGEN, 
and sorts them in a stack by increasing value. It then 
returns the stack and the number of elements in the stack. 
The flowchart and the listing are shown in Figs. C.10 
and C.ll, respectively. 



r r r o c, r, r. o r- 
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Reproduced from 
best available copy 



SUUKuUTlNE DOUBLE <L) 

THIS PROCEDURE SUBSTITUTES A PARTICULAR 2X2 MATRIX 
FOR EACH -l»Uf+l ELEMENT OF MATRIX H t THUS 
DOUBLING THE SIZE 

COMMON V(02>02) 

SET INDEXING SO THAT NO ORIGINAL ELEMENTS ARE ERASED 
DO 1 I - 1 » L 
II=L-I+1 
111=2*11 
DO 1 J= 1 » L 
00=L-0+l 
000—2*00 

TEST FOR A + I ELEMENT 
IMV(IIrOO)-O.Oi) if 2 '2 

2 V(III#000)=-i 
V(III»000-1)=1 
V ( I I I “ 1 r 000 ) = 1 
V( III-1»000-1)=1 
GO TO 1 

TEST FOR A -1 ELEMENT 

3 IE ( V ( 1 1 * 00 ) +0 • U I ) 4#4rS 

4 V(III»000)=1 
V(IIIfOOO-l)=-i 
V(IH-IfOOO)=-i 
V( III-1#000-1)=-1 
GO TO 1 

C 

0 ASSUME A 0 ELEMENT BY DEFAULT 

5 V ( 1 1 1 # 000 ) =-l 
V(IIi»OUO-l)=-l 
V(III-lrOOO)=-l 
V( III-l»000-*i)=l 

1 CONTINUE 

C 

RETURN 
END . 


Fig. C,9, Listing of Subroutine DOUBLE 



127 



Fig. C.10. Subroutine QADRES 
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SUljUuUT I NL GAUHtb < P • L * HI. *» ) 
e THIb PROCEDURE COMPUTES THE QUADRATIC RESIDUES 

^ OH ANY PRIME P AND STOKES THEM IN STACK RES(L) 

C 

INTEGER P » GRES » PM ♦ RES < 20 ) 

L-0 

PM=P-1 

DO 1 1=1 t PM 

1SG=I*I 

IF ( ISG • L T • P ) GO TO 1 

C COMPUTE TIC GUAURaTIC RESIDUE FOR EACH INTEGER 

GRES=ISQ-USG/P)*P 
IF (L.GT.O) GO TO 4 

rls(1)=qrf:s 

GO TO 3 

4 CONTINUE 

C CHECK THE STACK RES FOR DUPLICATION 

DO 5 J=1m_ 

IF (RES ( J) . KQ.GKES) GO TO 1 

5 CONTINUE 

L AUU THE NEW QUADRATIC RESIDUE TO THE STACK AND 

C SHIFT POINTER L 

RESIL+l )=GRES 
DO 2 J=I#L 

IF ( GRES . GT . RES ( L- J+l ) ) GO TO 3 
RtS(L~J+2)=RES(L-J+l) 

RLS(L-J+1)=GRES 
L=L + 1 
X CONTINUE 

RETURN 
END 


KP 


Fig. C.il Listing of Subroutine QADRES 
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C.3 COMPUTATION OF EIGENVALUES AND EIGENVECTORS 

The computer simulation described in this appendix 
uses one of two methods to obtain the eigenvalues and 
eigenvectors of an N x N real symmetric matrix A* Both 
methods are iterative, thus adaptable for machine compu- 
tation, and both make use of a series of unitary similarity 
transformations to reduce the given matrix to a simplified 
form. 

The first method [54] , originally proposed by Jacobi 
in 1864, has been revised in recent years for application 
to modern digital computers. The original method. is based 
on applying an infinite series of particular similarity 
transformations of the form A i+1 = with A Q = A. 

The matrices are plane (two-dimensional) rotations, and 
by properly selecting the plane and angle of rotation , any 
particular off-diagonal element can be annihilated, i.e., 
reduced to zero. Although the annihilation of other elements 
will in general cause previously annihilated elements to be 
nonzero, Jacobi showed that all off-diagonal elements tend 
to zero in the limit. By defining S_ and D to be the limits 
of S. and D. , respectively, it follows that 

D = S T A S , 

with D diagonal and S unitary. This implies that the 
diagonal elements of D are the eigenvalues of A, with 
the corresponding eigenvectors being the columns of S. 
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With machine computations the process is, of course, 
terminated after a finite number of steps as soon as the 
magnitude of the largest off-diagonal element becomes less 
than some predetermined threshold value. The convergence 
rate can also be improved by properly selecting the order 
in which the elements are annihilated. 

The second method 155] uses a procedure suggested 
by Householder to convert a real symmetric matrix into 
tridiagonal (codiagonal) form C with a series of N - 2 
unitary transformations of the form 1 - 2 u u T where | |H j [ =1 ^ 
By properly selecting u, any column and the corresponding 
row can be transformed to contain only two elements, thus 
giving the tridiagonal form. Once A is in tridiagonal 
form the eigenvalues of C, hence of A, are found by 
noticing that the set of characteristic polynomials 
corresponding to the leading principle submatrices forms 
a Sturm sequence. Each root of the largest polynomial is 
an eigenvalue of A, and the roots are found by obtaining 

successively smaller intervals in which each lies. The 
eigenvectors are then found for C by applying the inverse 
power method 156] which is particularly simple for 
codiagonal matrices. These eigenvectors are then transformed 
into those of A by using the same series of unitary trans- 
formations obtained in computing C initially. 
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